Advanced Statistics with R code 90% doen just check it and make a summary

profilearoojch
s_p500analysis_inr.pdf

S&P500 analysis Fitting an AR(1)-GARCH(1,1) model with Gaussian errors > # Fitting a ARMA(0,0)-GARCH(1,1) model using the rugarch package

> # log returns show weak serial autocorrelations in the mean and an ARCH

effect

> # Fitting an GARCH(1,1) model

> # Use ugarchspec() function to specify model

> garch11.spec=ugarchspec(variance.model=list(garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)))

> #estimate model

> garch11.fit=ugarchfit(spec=garch11.spec, data=ret)

> garch11.fit

*---------------------------------*

* GARCH Model Fit *

*---------------------------------*

Conditional Variance Dynamics

-----------------------------------

GARCH Model : sGARCH(1,1)

Mean Model : ARFIMA(0,0,0)

Distribution : norm

Optimal Parameters

------------------------------------

Estimate Std. Error t value Pr(>|t|)

mu 0.006814 0.001798 3.7889 0.000151

omega 0.000090 0.000049 1.8302 0.067224

alpha1 0.118390 0.031384 3.7723 0.000162

beta1 0.844396 0.034387 24.5556 0.000000

The fitted ARMA(0,0)-GARCH(1,1) model with Gaussian errors can be written as

222

8444.01184.000009.0

0068.0

tttttt

tt

aea

ar

 



Robust Standard Errors:

Estimate Std. Error t value Pr(>|t|)

mu 0.006814 0.001914 3.5592 0.000372

omega 0.000090 0.000070 1.2858 0.198502

alpha1 0.118390 0.041165 2.8760 0.004028

beta1 0.844396 0.040113 21.0503 0.000000

LogLikelihood : 828.5681

Information Criteria

------------------------------------

Akaike -3.4286

Bayes -3.3938

Shibata -3.4287

Hannan-Quinn -3.4149

Q-Statistics on Standardized Residuals

------------------------------------

statistic p-value

Lag[1] 0.5261 0.4682

Lag[p+q+1][1] 0.5261 0.4682

Lag[p+q+5][5] 5.9001 0.3161

d.o.f=0

H0 : No serial correlation

Q-Statistics on Standardized Squared Residuals

------------------------------------

statistic p-value

Lag[1] 0.2657 0.6062

Lag[p+q+1][3] 0.8252 0.3637

Lag[p+q+5][7] 3.0889 0.6863

d.o.f=2

ARCH LM Tests

------------------------------------

Statistic DoF P-Value

ARCH Lag[2] 0.5084 2 0.7755

ARCH Lag[5] 0.9231 5 0.9685

ARCH Lag[10] 6.3220 10 0.7875

Nyblom stability test

------------------------------------

Joint Statistic: 1.3254

Individual Statistics:

mu 0.07929

omega 0.20176

alpha1 0.02459

beta1 0.09300

Asymptotic Critical Values (10% 5% 1%)

Joint Statistic: 1.07 1.24 1.6

Individual Statistic: 0.35 0.47 0.75

Sign Bias Test

------------------------------------

t-value prob sig

Sign Bias 2.1382 0.0330120 **

Negative Sign Bias 0.5723 0.5673614

Positive Sign Bias 1.3403 0.1807809

Joint Effect 17.4223 0.0005786 ***

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 43.12 0.0012496

2 30 53.82 0.0033902

3 40 75.76 0.0003804

4 50 79.60 0.0037117

Ljung-Box test for serial correlation

computed on residuals, and Ljung-

Box test for ARCH/GARCH effect

computed on squared residuals.

Goodness of fit test for distribution of error

term. The null hypothesis states that the

distribution for the error terms in the model

is adequate. Thus a small p-value < α,

indicates that null hypothesis can be

rejected and the distribution assumption is

not adequate.

> #create selection list of plots for garch(1,1) fit

> plot(garch11.fit)

> #to display all subplots on one page

> plot(garch11.fit, which="all")

Residual analysis of GARCH model shows that the model fits the data adequately. Ljung Box test (Q-statistic) on residuals is not significant, showing that hypothesis of no correlation for residuals cannot be rejected. Similarly the Ljung Box test (Q-statistic) on the squared standardized residuals is not significant suggesting that residuals show no ARCH/GARCH effect. The Adjusted Pearson goodness of fit test is significant indicating that the normal distribution assumed for the error term is not appropriate, as also shown by the QQ plot of the residuals.

Fitting an AR(1)-GARCH(1,1) model with t-distributed errors

Elapsed time : 0.260026

> plot(garch11.fit, which="all")

Error in plot.new() : figure margins too large

> # use Student-t innovations

> #specify model using functions in rugarch package

> #Fit ARMA(0,0)-GARCH(1,1) model with t-distribution

> garch11.t.spec=ugarchspec(variance.model=list(garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)), distribution.model = "std")

> #estimate model

> garch11.t.fit=ugarchfit(spec=garch11.t.spec, data=ret)

> garch11.t.fit

*---------------------------------*

* GARCH Model Fit *

*---------------------------------*

Conditional Variance Dynamics

-----------------------------------

GARCH Model : sGARCH(1,1)

Mean Model : ARFIMA(0,0,0)

Distribution : std

Optimal Parameters

------------------------------------

Estimate Std. Error t value Pr(>|t|)

mu 0.008012 0.001746 4.5892 0.000004

omega 0.000122 0.000068 1.7901 0.073444

alpha1 0.123683 0.038919 3.1780 0.001483

beta1 0.821471 0.050706 16.2008 0.000000

shape 6.891814 1.978494 3.4834 0.000495

The fitted ARMA(0,0)-GARCH(1,1) model with Gaussian errors can be written as

222

8215.01237.000012.0

0080.0

tttttt

tt

aea

ar

 



Error term {et} has t-distribution with 7 degrees of freedom

Robust Standard Errors:

Estimate Std. Error t value Pr(>|t|)

mu 0.008012 0.001817 4.4090 0.000010

omega 0.000122 0.000071 1.7216 0.085136

alpha1 0.123683 0.038112 3.2453 0.001173

beta1 0.821471 0.041476 19.8061 0.000000

shape 6.891814 1.976630 3.4866 0.000489

LogLikelihood : 839.2684

Information Criteria

------------------------------------

Akaike -3.4689

Bayes -3.4255

Shibata -3.4691

Hannan-Quinn -3.4518

Q-Statistics on Standardized Residuals

------------------------------------

statistic p-value

Lag[1] 0.4662 0.4947

Lag[p+q+1][1] 0.4662 0.4947

Lag[p+q+5][5] 5.6618 0.3405

d.o.f=0

H0 : No serial correlation

Q-Statistics on Standardized Squared Residuals

------------------------------------

statistic p-value

Lag[1] 0.2474 0.6189

Lag[p+q+1][3] 0.8569 0.3546

Lag[p+q+5][7] 2.9570 0.7066

d.o.f=2

ARCH LM Tests

------------------------------------

Statistic DoF P-Value

ARCH Lag[2] 0.5959 2 0.7423

ARCH Lag[5] 0.9740 5 0.9646

ARCH Lag[10] 6.0871 10 0.8079

Nyblom stability test

------------------------------------

Joint Statistic: 1.396

Individual Statistics:

mu 0.17184

omega 0.23451

alpha1 0.04303

beta1 0.10209

shape 0.07244

Asymptotic Critical Values (10% 5% 1%)

Joint Statistic: 1.28 1.47 1.88

Individual Statistic: 0.35 0.47 0.75

Sign Bias Test

------------------------------------

t-value prob sig

Sign Bias 2.121 0.0343991 **

Negative Sign Bias 0.549 0.5832924

Positive Sign Bias 1.206 0.2284569

Joint Effect 16.352 0.0009603 ***

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 22.66 0.25269

2 30 44.47 0.03312

3 40 39.50 0.44759

4 50 43.64 0.68967

Residual analysis of GARCH model with t-distributed error terms shows that the model fits the data adequately. Ljung Box test (Q-statistic) on residuals is not significant, showing that hypothesis of no correlation for residuals cannot be rejected. Similarly the Ljung Box test (Q-statistic) on the squared standardized residuals is not significant suggesting that residuals show no ARCH/GARCH effect. The Adjusted Pearson goodness of fit test is not significant indicating that the distribution of the error terms can be described by a t-distribution with 7 degrees of freedom,

as also shown by the QQ plot of the residuals (created using plot(garch11.t.fit, which=9).

Fitting an ARMA(0,0)-EGARCH(1,1) model The EGARCH model fitted by the rugarch package has a slightly different form than the model in the textbook. Here is the general expression of the EGARCH(1,1) model:

)ln(|))(||(|()ln(

,

2

1111111

2

 



ttttt

tttttt

eEee

eaar





where µt can follow an ARMA process, but most typically it will be constant. Gamma1 (γ1) is the leverage parameter. So if gamma1 in output is significant, then we can conclude that the volatility process has an asymmetric behavior.

Fitting an ARMA(0,0)-EGARCH(1,1) model with Gaussian distribution (similar to SAS example)

> #Fit ARMA(0,0)-eGARCH(1,1) model with Gaussian distribution

> egarch11.spec=ugarchspec(variance.model=list(model = "eGARCH",

garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0)))

> #estimate model

> egarch11.fit=ugarchfit(spec=egarch11.spec, data=ret)

> egarch11.fit

*---------------------------------*

* GARCH Model Fit *

*---------------------------------*

Conditional Variance Dynamics

-----------------------------------

GARCH Model : eGARCH(1,1)

Mean Model : ARFIMA(0,0,0)

Distribution : norm

Optimal Parameters

------------------------------------

Estimate Std. Error t value Pr(>|t|)

mu 0.005937 0.001860 3.1915 0.001415

omega -0.610497 0.263901 -2.3134 0.020703

alpha1 -0.111069 0.042614 -2.6064 0.009150

beta1 0.902462 0.041695 21.6442 0.000000

gamma1 0.209405 0.050542 4.1432 0.000034

Using the R output above, the ARMA(0,0)-EGARCH(1,1) model can be written as follows.

Fitted model:

rt = 0.0059 + at, at=σtet

ln(σ2t) = -0.610 + (-0.111 et-1 + 0.209(|et-1| - E(|et-1|)) + 0.9024 ln(σ 2

t-1) Note that since et has Gaussian distribution, the E(|et|)=sqrt(2/pi) = 0.7979 or approx 0.80 (see page 143 in textbook). Thus we can rewrite the expression above as

  



 





07979.0209.0209.0111.0

07979.0209.0209.0111.0 )ln(9024.0610.0)ln(

111

1112

1

2

ttt

ttt

tt

eifee

eifee 

And after some algebra, we can write:

  



 





0320.0

0098.0 )ln(9024.07767.0)ln(

11

112

1

2

tt

tt

tt

eife

eife 

Taking the antilog transformation we have

  



 





0)320.0exp(

0)098.0exp( )7767.0exp(

11

119024.02

1

2

tt

ttx

tt

eife

eife 

For a standardized shock with magnitude 2, (i.e. two standard deviations), we have

56.1 )2098.0exp(

))2(320.0exp(

)2(

)2(

1

2

1

2

 

 



tt

tt

e

e

Therefore, the impact of negative shock of size two-standard deviations is about 56% higher than the impact of a positive shock of the same size.

Robust Standard Errors:

Estimate Std. Error t value Pr(>|t|)

mu 0.005937 0.001958 3.0328 0.002423

omega -0.610497 0.364615 -1.6744 0.094060

alpha1 -0.111069 0.066343 -1.6742 0.094096

beta1 0.902462 0.057336 15.7400 0.000000

gamma1 0.209405 0.048127 4.3511 0.000014

LogLikelihood : 835.9936

Information Criteria

------------------------------------

Akaike -3.4553

Bayes -3.4119

Shibata -3.4555

Hannan-Quinn -3.4382

Q-Statistics on Standardized Residuals

------------------------------------

statistic p-value

Lag[1] 0.3351 0.5627

Lag[p+q+1][1] 0.3351 0.5627

Lag[p+q+5][5] 5.1570 0.3970

d.o.f=0

H0 : No serial correlation

Q-Statistics on Standardized Squared Residuals

------------------------------------

statistic p-value

Lag[1] 0.7262 0.3941

Lag[p+q+1][3] 1.0203 0.3124

Lag[p+q+5][7] 2.8738 0.7194

d.o.f=2

ARCH LM Tests

------------------------------------

Statistic DoF P-Value

ARCH Lag[2] 1.018 2 0.6010

ARCH Lag[5] 1.165 5 0.9482

ARCH Lag[10] 6.026 10 0.8131

Nyblom stability test

------------------------------------

Joint Statistic: 1.5861

Individual Statistics:

mu 0.1829

omega 0.1224

alpha1 0.1527

beta1 0.1301

gamma1 1.0125

Asymptotic Critical Values (10% 5% 1%)

Joint Statistic: 1.28 1.47 1.88

Individual Statistic: 0.35 0.47 0.75

Sign Bias Test

------------------------------------

t-value prob sig

Sign Bias 1.973 0.04902 **

Negative Sign Bias 1.165 0.24478

Positive Sign Bias 1.041 0.29840

Joint Effect 11.259 0.01041 **

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 22.49 0.2604

2 30 35.86 0.1777

3 40 50.31 0.1060

4 50 55.07 0.2558

Residual analysis of EGARCH model shows that the model fits the data adequately – residuals are white noise and show no ARCH effect. The goodness of fit test supports the choice of a Gaussian distribution for the error term. Although the qqplot shows that the error distribution has thicker left tail than the normal distribution. We will fit the t-distribution to check if that’s a better fit for the behavior of extreme values.

Fitting an ARMA(0,0)-EGARCH(1,1) model with t-distribution

> #Fit ARMA(0,0)-eGARCH(1,1) model with t-distribution

> egarch11.t.spec=ugarchspec(variance.model=list(model = "eGARCH",

garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0)),

distribution.model = "std")

> #estimate model

> egarch11.t.fit=ugarchfit(spec=egarch11.t.spec, data=ret)

> egarch11.t.fit

*---------------------------------*

* GARCH Model Fit *

*---------------------------------*

Conditional Variance Dynamics

-----------------------------------

GARCH Model : eGARCH(1,1)

Mean Model : ARFIMA(0,0,0)

Distribution : std

Optimal Parameters

------------------------------------

Estimate Std. Error t value Pr(>|t|)

mu 0.007093 0.001757 4.0369 0.000054

omega -0.677075 0.246000 -2.7523 0.005917

alpha1 -0.145367 0.045752 -3.1773 0.001487

beta1 0.893975 0.038715 23.0912 0.000000

gamma1 0.202717 0.055934 3.6242 0.000290

shape 7.926664 2.405188 3.2957 0.000982

Using the R output above, the ARMA(0,0)-EGARCH(1,1) model can be written as follows.

Fitted model:

rt = 0.007 + at, at=σtet

ln(σ2t) = -0.677 + (-0.145 et-1 + 0.203(|et-1| - E(|et-1|)) + 0.894 ln(σ 2

t-1) with t-distribution with 8 degrees of freedom (nearest integer to shape value)

Note that since et has t-distribution, the E(|et|)= 



)2/()1(

]2/)1[(22



 (see page 143 in book) where ν=d.f. of t-

distribution (denoted by shape in R output). Note that since gamma1 is significant, the volatility has an asymmetric behavior. Therefore a negative shock has a stronger impact on the volatility compared to a positive shock of the same size.

Robust Standard Errors:

Estimate Std. Error t value Pr(>|t|)

mu 0.007093 0.001888 3.7574 0.000172

omega -0.677075 0.212770 -3.1822 0.001462

alpha1 -0.145367 0.042936 -3.3857 0.000710

beta1 0.893975 0.034112 26.2071 0.000000

gamma1 0.202717 0.045692 4.4366 0.000009

shape 7.926664 2.720128 2.9141 0.003567

LogLikelihood : 846.92

Information Criteria

------------------------------------

Akaike -3.4965

Bayes -3.4445

Shibata -3.4969

Hannan-Quinn -3.4761

Q-Statistics on Standardized Residuals

------------------------------------

statistic p-value

Lag[1] 0.1872 0.6653

Lag[p+q+1][1] 0.1872 0.6653

Lag[p+q+5][5] 4.7776 0.4436

d.o.f=0

H0 : No serial correlation

Q-Statistics on Standardized Squared Residuals

------------------------------------

statistic p-value

Lag[1] 0.6769 0.4107

Lag[p+q+1][3] 0.9526 0.3291

Lag[p+q+5][7] 2.5747 0.7652

d.o.f=2

ARCH LM Tests

------------------------------------

Statistic DoF P-Value

ARCH Lag[2] 0.9725 2 0.6149

ARCH Lag[5] 1.0934 5 0.9547

ARCH Lag[10] 5.4934 10 0.8559

Nyblom stability test

------------------------------------

Joint Statistic: 1.6436

Individual Statistics:

mu 0.20752

omega 0.12081

alpha1 0.15261

beta1 0.13304

gamma1 0.93741

shape 0.08736

Asymptotic Critical Values (10% 5% 1%)

Joint Statistic: 1.49 1.68 2.12

Individual Statistic: 0.35 0.47 0.75

Sign Bias Test

------------------------------------

t-value prob sig

Sign Bias 1.7639 0.07840 *

Negative Sign Bias 1.2413 0.21510

Positive Sign Bias 0.9747 0.33022

Joint Effect 8.8291 0.03165 **

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 20.66 0.35571

2 30 39.23 0.09736

3 40 54.47 0.05098

4 50 66.51 0.04860

Residual analysis of EGARCH model shows that the model fits the data adequately – residuals are white noise and show no ARCH effect. However the goodness of fit test shows that the t-distribution is not a good choice for the error terms (test p-values are small and rejects the null hypothesis of error term having a t-distribution. The qqplot shows that the error distribution has thicker tails than the t-distribution. (Further analysis shows that the ged distribution does a better job representing the extreme values distribution)

Fitting an ARMA(0,0)-GJRGARCH(1,1) model or TGARCH model

The GJRGARCH(1,1) model fitted by the rugarch package has the following expression

)ln()()ln(

,

2

11

2

1111

2

 



tttt

tttttt

aN

eaar





where µt can follow an ARMA process, but most typically it will be constant. Nt-1 is the indicator variable s.t. Nt-1= 1 when at-1 (shock at time t-1) is negative, and Nt-1 = 0 otherwise. Gamma1 (γ1) is the leverage parameter. So if gamma1 in output is significant, then we can conclude that the volatility process has an asymmetric behavior. > #Fit ARMA(0,0)-TGARCH(1,1) model with t-distribution

> gjrgarch11.t.spec=ugarchspec(variance.model=list(model = "gjrGARCH",

garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0)), distribution.model = "std")

> #estimate model

> gjrgarch11.t.fit=ugarchfit(spec=gjrgarch11.t.spec, data=ret)

Warning message:

In .makefitmodel(garchmodel = "gjrGARCH", f = .gjrgarchLLH, T = T, :

NaNs produced

> gjrgarch11.t.fit

*---------------------------------*

* GARCH Model Fit *

*---------------------------------*

Conditional Variance Dynamics

-----------------------------------

GARCH Model : gjrGARCH(1,1)

Mean Model : ARFIMA(0,0,0)

Distribution : std

Optimal Parameters

------------------------------------

Estimate Std. Error t value Pr(>|t|)

mu 0.007202 0.001766 4.077739 0.000045

omega 0.000212 0.000094 2.262826 0.023646

alpha1 0.000001 0.001394 0.000791 0.999369

beta1 0.781111 0.065122 11.994488 0.000000

gamma1 0.215386 0.069611 3.094141 0.001974

shape 7.179979 2.031526 3.534279 0.000409

Using the R output above, the ARMA(0,0)-TGARCH(1,1) model can be written as follows. Note that the arch(1) coefficient is zero and we remove it from the model.

2

1

2

11

2 781.0)215.00000.0(00021.0

,0072.0

 



tttt

ttttt

aN

eaar



Where error term is assumed to have t-distribution with 7 degrees of freedom.

or

2

1

2

11

2 781.0215.000021.0

,0072.0

 



tttt

ttttt

aN

eaar



Where Nt-1 is an indicator variable for negative innovations such that   

 

00

01

1

1

1

t

t

t

aif

aif N

For a standardized shock with magnitude 2, (i.e. two standard deviations), we have

966.1 002.0781.0)002.04()000.0(00021.0

002.0781.0)002.04()215.00000.0(00021.0

)2(

)2(

1

2

1

2

 

 



tt

tt

e

e

Where the value for 2

1t a is computed as

2

1

22

1  

ttt ea  . Since

2

1t  in the expression above is unknown, we set

2

1t  = sample variance of the S&P500 index returns (that is

2

1t  = 0.002) .

Robust Standard Errors:

Estimate Std. Error t value Pr(>|t|)

mu 0.007202 0.001921 3.749767 0.000177

omega 0.000212 0.000112 1.902791 0.057068

alpha1 0.000001 0.000054 0.020238 0.983853

beta1 0.781111 0.063980 12.208691 0.000000

gamma1 0.215386 0.056979 3.780097 0.000157

shape 7.179979 2.274951 3.156102 0.001599

LogLikelihood : 844.6336

Information Criteria

------------------------------------

Akaike -3.4870

Bayes -3.4350

Shibata -3.4873

Hannan-Quinn -3.4666

Q-Statistics on Standardized Residuals

------------------------------------

statistic p-value

Lag[1] 0.1348 0.7135

Lag[p+q+1][1] 0.1348 0.7135

Lag[p+q+5][5] 4.7791 0.4434

d.o.f=0

H0 : No serial correlation

Q-Statistics on Standardized Squared Residuals

------------------------------------

statistic p-value

Lag[1] 0.7412 0.3893

Lag[p+q+1][3] 1.0145 0.3138

Lag[p+q+5][7] 2.5887 0.7631

d.o.f=2

ARCH LM Tests

------------------------------------

Statistic DoF P-Value

ARCH Lag[2] 1.015 2 0.6020

ARCH Lag[5] 1.239 5 0.9411

ARCH Lag[10] 5.120 10 0.8830

Nyblom stability test

------------------------------------

Joint Statistic: 2.2942

Individual Statistics:

mu 0.18155

omega 0.27666

alpha1 0.05641

beta1 0.20413

gamma1 0.03198

shape 0.06486

Asymptotic Critical Values (10% 5% 1%)

Joint Statistic: 1.49 1.68 2.12

Individual Statistic: 0.35 0.47 0.75

Sign Bias Test

------------------------------------

t-value prob sig

Sign Bias 2.117 0.03474 **

Negative Sign Bias 1.408 0.15987

Positive Sign Bias 0.766 0.44408

Joint Effect 10.173 0.01715 **

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 19.83 0.4048

2 30 29.37 0.4457

3 40 37.17 0.5535

4 50 52.78 0.3300

Elapsed time : 0.458046

> plot(gjrgarch11.t.fit, which="all")

Residual analysis of TGARCH model shows that the model fits the data adequately – residuals are white noise and show no ARCH effect. The goodness of fit test also shows that the t-distribution is adequate.

Apply information criteria for model selection

> # MODEL COMPARISON

> # compare information criteria

> model.list = list(garch11 = garch11.fit, garch11.t = garch11.t.fit,

+ egarch11 = egarch11.t.fit,

+ gjrgarch11 = gjrgarch11.t.fit)

> info.mat = sapply(model.list, infocriteria)

> rownames(info.mat) = rownames(infocriteria(garch11.fit))

> info.mat

garch11 garch11.t egarch11 gjrgarch11

Akaike -3.428558 -3.468892 -3.480644 -3.487042

Bayes -3.393831 -3.425483 -3.428554 -3.434952

Shibata -3.428694 -3.469105 -3.480950 -3.487348

Hannan-Quinn -3.414909 -3.451830 -3.460170 -3.466568

Best model according to model selection criteria is the GJR-GARCH(1,1) model with t-distribution R CODE: # Analysis of daily S&P500 index

#

library(fBasics)

library(tseries)

library(rugarch)

# import data in R

# import libraries for TS analysis

myd= read.table('sp500_feb1970_Feb2010.txt', header=T)

# create time series object •

rts= ts(myd$return, start = c(1970, 1), frequency=12)

# create a simple numeric object

ret =myd$return;

# CREATE TIME PLOT

plot(rts)

# Plots ACF function of vector data

acf(ret)

# Plot ACF of squared data to check for non-linear dependence

acf(ret^2)

# Computes Ljung-Box test on squared returns to test non-linear independence at lag 6

and 12

Box.test(ret^2,lag=6,type='Ljung')

Box.test(ret^2,lag=12,type='Ljung')

# Computes Ljung-Box test on absolute returns to test non-linear independence at lag 6

and 12

Box.test(abs(ret),lag=6,type='Ljung')

Box.test(abs(ret),lag=12,type='Ljung')

# FITTING AN GARCH(1,1) MODEL WITH GAUSSIAN DISTRIBUTION

# Use ugarchspec() function to specify model

garch11.spec=ugarchspec(variance.model=list(garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)))

#estimate model

garch11.fit=ugarchfit(spec=garch11.spec, data=ret)

garch11.fit

#persistence = alpha1+beta1

persistence(garch11.fit)

#half-life: ln(0.5)/ln(alpha1+beta1)

halflife(garch11.fit)

#create selection list of plots for garch(1,1) fit

plot(garch11.fit)

#to display all subplots on one page

plot(garch11.fit, which="all")

#FIT ARMA(0,0)-GARCH(1,1) MODEL WITH T-DISTRIBUTION

# specify model using functions in rugarch package

garch11.t.spec=ugarchspec(variance.model=list(garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)), distribution.model = "std")

#estimate model

garch11.t.fit=ugarchfit(spec=garch11.t.spec, data=ret)

garch11.t.fit

plot(garch11.t.fit)

#FIT ARMA(0,0)-EGARCH(1,1) MODEL WITH GAUSSIAN DISTRIBUTION

egarch11.spec=ugarchspec(variance.model=list(model = "eGARCH", garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)))

#estimate model

egarch11.fit=ugarchfit(spec=egarch11.spec, data=ret)

egarch11.fit

plot(egarch11.fit, which="all")

#FIT ARMA(0,0)-EGARCH(1,1) MODEL WITH T-DISTRIBUTION

egarch11.t.spec=ugarchspec(variance.model=list(model = "eGARCH", garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)), distribution.model = "ged")

#estimate model

egarch11.t.fit=ugarchfit(spec=egarch11.t.spec, data=ret)

egarch11.t.fit

plot(egarch11.t.fit, which="all")

# compute expected value E(|e|)

shape=coef(egarch11.t.fit)[6]

exp.abse=(2*sqrt(shape-2)*gamma((shape+1)/2))/((shape-1)*gamma(shape/2)*sqrt(pi))

#FIT ARMA(0,0)-TGARCH(1,1) MODEL WITH T-DISTRIBUTION

gjrgarch11.t.spec=ugarchspec(variance.model=list(model = "gjrGARCH",

garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0)), distribution.model = "std")

#estimate model

gjrgarch11.t.fit=ugarchfit(spec=gjrgarch11.t.spec, data=ret)

gjrgarch11.t.fit

plot(gjrgarch11.t.fit, which="all")

#FIT ARMA(0,0)-IGARCH(1,1) MODEL WITH SKEWED T-DISTRIBUTION

igarch11.t.spec=ugarchspec(variance.model=list(model = "iGARCH", garchOrder=c(1,1)),

mean.model=list(armaOrder=c(0,0)), distribution.model = "std")

#estimate model

igarch11.t.fit=ugarchfit(spec=igarch11.t.spec, data=ret)

igarch11.t.fit

plot(igarch11.t.fit, which="all")

# MODEL COMPARISON

# compare information criteria

model.list = list(garch11 = garch11.fit, garch11.t = garch11.t.fit,

egarch11 = egarch11.t.fit,

gjrgarch11 = gjrgarch11.t.fit)

info.mat = sapply(model.list, infocriteria)

rownames(info.mat) = rownames(infocriteria(garch11.fit))

info.mat

# RE-FIT MODELS LEAVING 100 OUT-OF-SAMPLE OBSERVATIONS FOR FORECAST

# EVALUATION STATISTICS

garch11.fit = ugarchfit(spec=garch11.spec, data=ret, out.sample=100)

garch11.t.fit = ugarchfit(spec=garch11.t.spec, data=ret, out.sample=100)

egarch11.t.fit = ugarchfit(egarch11.t.spec, data=ret, out.sample=100)

tgarch11.t.fit = ugarchfit(spec=gjrgarch11.t.spec, data=ret, out.sample=100)

# COMPUTE 100 1-STEP AHEAD ROLLING FORECASTS W/O RE-ESTIMATING

garch11.fcst = ugarchforecast(garch11.fit, n.roll=100, n.ahead=1)

garch11.t.fcst = ugarchforecast(garch11.t.fit, n.roll=100, n.ahead=1)

egarch11.t.fcst = ugarchforecast(egarch11.t.fit, n.roll=100, n.ahead=1)

tgarch11.t.fcst = ugarchforecast(tgarch11.t.fit, n.roll=100, n.ahead=1)

# COMPUTE FORECAST EVALUATION STATISTICS USING FPM() FUNCTION

fcst.list = list(garch11=garch11.fcst, garch11.t=garch11.t.fcst,

egarch11.t=egarch11.t.fcst,

tgarch11.t=tgarch11.t.fcst)

fpm.mat = sapply(fcst.list, fpm)

fpm.mat