Advanced Statistics with R code 90% doen just check it and make a summary
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
1t a is computed as
2
1
22
1
ttt ea . Since
2
1t in the expression above is unknown, we set
2
1t = sample variance of the S&P500 index returns (that is
2
1t = 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