Homework 4

Problem 2

The first step is to import Microsoft stock price, GOOG/NASDAQ_MSFT, using Quandl packages. We will denote GOOG/NASDAQ_MSF as \(MSFT_{t}\).

MSFT <- Quandl("GOOG/NASDAQ_MSFT", type="zoo")
msft <- MSFT[,4]

Next step is to construct the log return of Microsoft stock, \(y_{t}\), as follow:

yt <- diff(log(msft))

The following figures shows the original time series of MSFT closing price and the log return, \(y_{t}\).

par(mfrow=c(1,2), cex=0.5, mar=c(2,2,3,1))
plot(msft, xlab="", ylab="MSFT Closing Price", main="Microsoft Stock Closing Price")
plot(yt, xlab="", ylab="", main="Log Return of MSFT")

(a) ARMA Model Estimation

The first step in model identification is to plot the ACF and PACF for all transformation of \(y_{t}\):

From the ACF and PACF of \(y_{t}\), the most likely specification is ARIMA(4,0,2), next, We can estimate the specification using the following code:

m1 <- arima(coredata(yt),order=c(4,0,5))
m1

## 
## Call:
## arima(x = coredata(yt), order = c(4, 0, 5))
## 
## Coefficients:
##           ar1     ar2     ar3     ar4     ma1      ma2      ma3     ma4
##       -0.3935  0.0001  0.1008  -0.601  0.3553  -0.0617  -0.1535  0.5802
## s.e.   0.1833  0.1603  0.1675   0.175  0.1835   0.1577   0.1727  0.1794
##           ma5  intercept
##       -0.0284      8e-04
## s.e.   0.0183      2e-04
## 
## sigma^2 estimated as 0.0005343:  log likelihood = 17810.1,  aic = -35598.2

tsdiag(m1,gof.lag=12)

Based on the plots above, ARIMA(4,0,5)specification is adequate. We can also test the specification further for adequacy using Ljung-Box test on the residuals. Given that we have 7582 observations, \(\log{T}\) is 8.9335 or approximately 8, we choose \(m=8\).

Box.test(coredata(m1$residuals), lag=8, type="Ljung")

## 
##  Box-Ljung test
## 
## data:  coredata(m1$residuals)
## X-squared = 3.0272, df = 8, p-value = 0.9326

From the estimation result above, we can confirm that ARIMA(4,0,5) is indeed adequate.

(b) ARCH(m) Model Estimation with Normal Innovations

Testing for ARCH Effect

First step is to specify a mean equation and test it for serial dependence, we define a new series dyt as squared demeaned form of \(y_{t}\):

dyt <- (coredata(yt)-mean(yt))^2
plot(coredata(dyt), type="l", main="MSFT squared demeaned daily log return")

Next we test for serial dependence for dyt as follow:

Box.test(coredata(yt), lag=8, type="Ljung")

## 
##  Box-Ljung test
## 
## data:  coredata(yt)
## X-squared = 41.78, df = 8, p-value = 1.49e-06

The \(Q(m)\) statistics of \(y_{t}\) is \(Q(8)=41.774\) with a \(p\)-value that close to zero, which mean that we have a serial dependence, hence we cannot test for ARCH effect using \((y_{t} - \bar{y_{t}})^2\) transformation.

Alternatively, we can use the square residual of the previously estimated ARMA model of \(y_{t}\) in part (a), i.e. we create a new series \((\hat{y_{t}} - \bar{y_{t}})^2\). One advantage to test post estimation residual of the ARMA model is that we already removed any linear dependence in the series.

Box.test(coredata(m1$residuals^2), lag=8, type="Ljung")

## 
##  Box-Ljung test
## 
## data:  coredata(m1$residuals^2)
## X-squared = 671.05, df = 8, p-value < 2.2e-16

From the test above, Ljung-Box statistics of the \((\hat{y_{t}} - \bar{y_{t}})^2\) shows strong ARCH effects with \(Q(8)=670.08\), the \(p\)-value of which is close to zero. Now we can proceed to the next step of ARCH order determination

Order Determination

In this part we will use the ACF and PACF of the series \((\hat{y_{t}} - \bar{y_{t}})^2\) to identify the order of ARCH.

From the PACF of \((\hat{y_{t}} - \bar{y_{t}})^2\) above, we can see that we have some big spikes that are stastically significant for lag 1 until lag 8. This indicate that we have ARCH(8) identified from the squared residuals of the ARMA(4,0,5).

Hence, now we can estimate the ARCH model using the following code:

m2 <- garchFit(~ arma(4,5) + garch(8,0), data = coredata(yt), trace = FALSE)
summary(m2)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 5) + garch(8, 0), data = coredata(yt), 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 5) + garch(8, 0)
## <environment: 0x00000000198ff178>
##  [data = coredata(yt)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  0.00025464   0.28588086   0.99129512   0.24176170  -0.79359851  
##         ma1          ma2          ma3          ma4          ma5  
## -0.31183399  -0.99999999  -0.25384343   0.81334942   0.00762050  
##       omega       alpha1       alpha2       alpha3       alpha4  
##  0.00017892   0.16032407   0.12407071   0.06421123   0.11464131  
##      alpha5       alpha6       alpha7       alpha8  
##  0.05573980   0.03014863   0.08959319   0.05416470  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      2.546e-04   6.142e-05    4.146 3.38e-05 ***
## ar1     2.859e-01          NA       NA       NA    
## ar2     9.913e-01          NA       NA       NA    
## ar3     2.418e-01   6.215e-02    3.890   0.0001 ***
## ar4    -7.936e-01   5.354e-02  -14.822  < 2e-16 ***
## ma1    -3.118e-01          NA       NA       NA    
## ma2    -1.000e+00          NA       NA       NA    
## ma3    -2.538e-01   6.185e-02   -4.104 4.06e-05 ***
## ma4     8.133e-01   5.905e-02   13.774  < 2e-16 ***
## ma5     7.620e-03   1.695e-02    0.450   0.6530    
## omega   1.789e-04   7.059e-06   25.345  < 2e-16 ***
## alpha1  1.603e-01   1.507e-02   10.639  < 2e-16 ***
## alpha2  1.241e-01   1.620e-02    7.658 1.89e-14 ***
## alpha3  6.421e-02   1.241e-02    5.173 2.30e-07 ***
## alpha4  1.146e-01   1.449e-02    7.914 2.44e-15 ***
## alpha5  5.574e-02   1.153e-02    4.834 1.34e-06 ***
## alpha6  3.015e-02   1.178e-02    2.559   0.0105 *  
## alpha7  8.959e-02   1.431e-02    6.262 3.81e-10 ***
## alpha8  5.416e-02   1.182e-02    4.584 4.56e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18649.88    normalized:  2.459108 
## 
## Description:
##  Thu Mar 24 12:56:47 2016 by user: Rullan Rinaldi 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  13095.76  0           
##  Shapiro-Wilk Test  R    W      NA        NA          
##  Ljung-Box Test     R    Q(10)  5.971373  0.8176626   
##  Ljung-Box Test     R    Q(15)  8.180336  0.9163493   
##  Ljung-Box Test     R    Q(20)  11.15978  0.9419731   
##  Ljung-Box Test     R^2  Q(10)  6.797968  0.7443704   
##  Ljung-Box Test     R^2  Q(15)  26.0947   0.037039    
##  Ljung-Box Test     R^2  Q(20)  53.63037  6.561534e-05
##  LM Arch Test       R    TR^2   7.02721   0.8558102   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -4.913206 -4.895835 -4.913218 -4.907244

(c) Model Checking for thE ARCH(8) Model

In this part we will check the adequacy of the previously estimated ARCH(8) model under the premises that if the model is properly specified, standardized residuals should be iid, thus they

• should not be serially correlated
• should not exhibit periods of increased volatility
• its squares standardized residuals should not be serially correlated

The first step is to extract the residuals from ARCH(8) estimation and transform it into standardized and squared standardized residuals as follow:

par(mfrow=c(2,2), cex=0.6, mar=c(2,2,3,1))
plot(m2@residuals, type="l", xlab="", ylab="", main="residuals")
plot(m2@sigma.t, type="l", xlab="", ylab="", main="conditional standard deviation")

stdresid <- m2@residuals/m2@sigma.t
sqrstdresid <- (m2@residuals/m2@sigma.t)^2

plot(stdresid, type="l", xlab="", ylab="", main="standardized residuals")
plot(sqrstdresid, type="l", xlab="", ylab="", main="squared standardized residuals")

The next step is to plot ACF and PACF for the standardized and squared standardized residuals as follow:

par(mfcol=c(2,2), cex=0.6, mar=c(2,2,3,1))
acf(stdresid, type="correlation", lag.max=24, ylab="", main="ACF for standardized residuals")
acf(stdresid, type="partial", lag.max=24, ylab="", main="PACF for standardized residuals")
acf(sqrstdresid, type="correlation", lag.max=24, ylab="", main="ACF for squared standardized residuals")
acf(sqrstdresid, type="partial", lag.max=24, ylab="", main="PACF for squared standardized residuals")

From the ACF for both standardized and squared standardized residuals, there is no spike at any lag except for lag 0 (which is always equal to one), this indicate that there is no serial correlation in either series.

To confirm this finding, we can conduct Ljung-Box test at selected lag (10, 15 and 20), to obtain the \(Q(10)\) and \(Q(20)\), we will use the following code for squared standardized residuals:

Box.test(sqrstdresid, lag = 10, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sqrstdresid
## X-squared = 6.798, df = 10, p-value = 0.7444

Box.test(sqrstdresid, lag = 15, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sqrstdresid
## X-squared = 26.095, df = 15, p-value = 0.03704

Box.test(sqrstdresid, lag = 20, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sqrstdresid
## X-squared = 53.63, df = 20, p-value = 6.562e-05

Jarque-Berra test with p-value of exactly zero indicate that the standardized residuals is normally distributed. From the Ljung-Box test we can see that at \(Q(10)\), the model is actually adequate, but caution have to be exercised, since at further lag (15 and 20) the result of Ljung-Box test is not quite good. This indicate that ARCH(8) probably not the best model, forecast based on this model is not suggested.

plot(m2, which=13)

The QQ-plot indicate that the residuals from our previously estimated model is actually normally distributed, however caution have to exercised since the tail of the distribution is deviating from the center line.

(d) GARCH(m,s) Model Estimation with Student-t Innovations

We will start our estimation with the most commonly used GARCH model with m=1 and s=1, i.e. GARCH(1,1) model. The estimation implemented in R using the following code:

m3 <- garchFit(~ arma(4,5) + garch(1,1), data = coredata(yt), trace = FALSE, cond.dist="sstd")
summary(m3)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 5) + garch(1, 1), data = coredata(yt), 
##     cond.dist = "sstd", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 5) + garch(1, 1)
## <environment: 0x000000001cb00db8>
##  [data = coredata(yt)]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  5.0814e-04   2.0009e-01   7.5285e-01   3.0909e-01  -9.2262e-01  
##         ma1          ma2          ma3          ma4          ma5  
## -2.3359e-01  -7.5191e-01  -2.9419e-01   9.3031e-01  -3.1521e-02  
##       omega       alpha1        beta1         skew        shape  
##  2.6024e-06   6.5364e-02   9.3366e-01   1.0410e+00   4.8014e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      5.081e-04   1.214e-04    4.187 2.83e-05 ***
## ar1     2.001e-01   3.957e-03   50.562  < 2e-16 ***
## ar2     7.528e-01   1.117e-02   67.427  < 2e-16 ***
## ar3     3.091e-01   1.339e-02   23.091  < 2e-16 ***
## ar4    -9.226e-01   9.451e-03  -97.623  < 2e-16 ***
## ma1    -2.336e-01   1.122e-02  -20.817  < 2e-16 ***
## ma2    -7.519e-01   9.484e-03  -79.284  < 2e-16 ***
## ma3    -2.942e-01   5.332e-03  -55.171  < 2e-16 ***
## ma4     9.303e-01   4.363e-03  213.227  < 2e-16 ***
## ma5    -3.152e-02   9.529e-03   -3.308 0.000940 ***
## omega   2.602e-06   7.317e-07    3.557 0.000376 ***
## alpha1  6.536e-02   8.764e-03    7.458 8.79e-14 ***
## beta1   9.337e-01   8.438e-03  110.653  < 2e-16 ***
## skew    1.041e+00   1.498e-02   69.510  < 2e-16 ***
## shape   4.801e+00   2.662e-01   18.033  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  19283.73    normalized:  2.542686 
## 
## Description:
##  Fri Mar 25 11:02:15 2016 by user: Rullan Rinaldi 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  24648.55  0        
##  Shapiro-Wilk Test  R    W      NA        NA       
##  Ljung-Box Test     R    Q(10)  5.02122   0.889756 
##  Ljung-Box Test     R    Q(15)  7.491683  0.9425441
##  Ljung-Box Test     R    Q(20)  10.16242  0.9651303
##  Ljung-Box Test     R^2  Q(10)  1.810779  0.9975954
##  Ljung-Box Test     R^2  Q(15)  7.094706  0.9549682
##  Ljung-Box Test     R^2  Q(20)  8.379393  0.9890451
##  LM Arch Test       R    TR^2   2.840243  0.9965706
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -5.081416 -5.067702 -5.081424 -5.076710

Next we extract the residuals from the estimation and transform it into standardized and squared standardized residuals as follow

stdresid2 <- m3@residuals/m2@sigma.t
sqrstdresid2 <- (m3@residuals/m2@sigma.t)^2
par(mfcol=c(2,2), cex=0.6, mar=c(2,2,3,1))
acf(stdresid2, type="correlation", lag.max=24, ylab="", main="ACF for standardized residuals")
acf(stdresid2, type="partial", lag.max=24, ylab="", main="PACF for standardized residuals")
acf(sqrstdresid2, type="correlation", lag.max=24, ylab="", main="ACF for squared standardized residuals")
acf(sqrstdresid2, type="partial", lag.max=24, ylab="", main="PACF for squared standardized residuals")

Plotting the QQ plot:

plot(m3, which=13)

The QQ-plot indicate that the residuals from our previously estimated model is actually normally distributed, however caution have to exercised since the tail of the distribution is deviating from the center line.

We also will estimate other frequently used specification of GARCH, the result can be seen in later part using stargazer packages. We implement the estimation using the following code:

m4 <- garchFit(~ arma(4,5) + garch(1,2), data = coredata(yt), trace = FALSE, cond.dist="sstd")
m5 <- garchFit(~ arma(4,5) + garch(2,1), data = coredata(yt), trace = FALSE, cond.dist="sstd")
m6 <- garchFit(~ arma(4,5) + garch(2,2), data = coredata(yt), trace = FALSE, cond.dist="sstd")

(e) TGARCH(m,s) Model Estimation

We will start our estimation with the most commonly used TGARCH model with m=1 and s=1, i.e. TGARCH(1,1) model. The estimation implemented in R using the following code:

par(mfrow=c(2,2))
m7 <- garchFit(~ arma(4,5) + garch(1,1), data=coredata(yt), trace=FALSE, leverage=TRUE)
summary(m7)

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 5) + garch(1, 1), data = coredata(yt), 
##     leverage = TRUE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 5) + garch(1, 1)
## <environment: 0x0000000022863ed8>
##  [data = coredata(yt)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  1.6727e-04   3.2114e-01   8.7109e-01   4.5649e-01  -9.0305e-01  
##         ma1          ma2          ma3          ma4          ma5  
## -3.4121e-01  -8.9277e-01  -4.5230e-01   9.2462e-01   2.9746e-03  
##       omega       alpha1       gamma1        beta1  
##  5.1112e-06   6.0272e-02   1.0623e-01   9.3098e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      1.673e-04   5.177e-05    3.231 0.001233 ** 
## ar1     3.211e-01          NA       NA       NA    
## ar2     8.711e-01          NA       NA       NA    
## ar3     4.565e-01   1.372e-02   33.266  < 2e-16 ***
## ar4    -9.030e-01   9.174e-03  -98.437  < 2e-16 ***
## ma1    -3.412e-01   7.673e-03  -44.468  < 2e-16 ***
## ma2    -8.928e-01          NA       NA       NA    
## ma3    -4.523e-01   8.478e-03  -53.350  < 2e-16 ***
## ma4     9.246e-01   2.030e-03  455.392  < 2e-16 ***
## ma5     2.975e-03   4.730e-03    0.629 0.529447    
## omega   5.111e-06   7.682e-07    6.654 2.86e-11 ***
## alpha1  6.027e-02   5.557e-03   10.846  < 2e-16 ***
## gamma1  1.062e-01   3.063e-02    3.468 0.000524 ***
## beta1   9.310e-01   6.127e-03  151.939  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18840.55    normalized:  2.484249 
## 
## Description:
##  Thu Mar 24 12:58:30 2016 by user: Rullan Rinaldi 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  13902     0        
##  Shapiro-Wilk Test  R    W      NA        NA       
##  Ljung-Box Test     R    Q(10)  2.360554  0.9927557
##  Ljung-Box Test     R    Q(15)  4.44773   0.9958595
##  Ljung-Box Test     R    Q(20)  7.851727  0.9928041
##  Ljung-Box Test     R^2  Q(10)  2.059226  0.9958657
##  Ljung-Box Test     R^2  Q(15)  7.186648  0.9522618
##  Ljung-Box Test     R^2  Q(20)  8.642609  0.9866892
##  LM Arch Test       R    TR^2   3.002204  0.9955284
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -4.964806 -4.952007 -4.964813 -4.960413

stdresid3 <- m7@residuals/m2@sigma.t
sqrstdresid3 <- (m7@residuals/m2@sigma.t)^2
par(mfcol=c(2,2), cex=0.6, mar=c(2,2,3,1))
acf(stdresid3, type="correlation", lag.max=24, ylab="", main="ACF for standardized residuals")
acf(stdresid3, type="partial", lag.max=24, ylab="", main="PACF for standardized residuals")
acf(sqrstdresid3, type="correlation", lag.max=24, ylab="", main="ACF for squared standardized residuals")
acf(sqrstdresid3, type="partial", lag.max=24, ylab="", main="PACF for squared standardized residuals")

Plotting the QQ plot:

plot(m7, which=13)

The QQ-plot indicate that the residuals from our previously estimated model is actually normally distributed, however caution have to exercised since the tail of the distribution is deviating from the center line.

We also will estimate other frequently used specification of GARCH, the result can be seen in later part using stargazer packages. We implement the estimation using the following code:

m8 <- garchFit(~ arma(4,5) + garch(1,2), data = coredata(yt), trace = FALSE, leverage=TRUE)
m9 <- garchFit(~ arma(4,5) + garch(2,1), data = coredata(yt), trace = FALSE, leverage=TRUE)
m10 <- garchFit(~ arma(4,5) + garch(2,2), data = coredata(yt), trace = FALSE, leverage=TRUE)

(e) Comparing the Results using StargazerX

First we need to set the working directory and load the function manually:

setwd("C:/Users/Rullan Rinaldi/Google Drive/ECO - 5316/HW 4")
source("stargazerX.r")
source("stargazerX-internal.r")

Constructing out comparison table using StargazerX

stargazerX(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10, type="html", title="Estimation Results for Log Return of MSFT", 
           align=TRUE, header=FALSE, model.numbers=TRUE,
           dep.var.caption="", 
           column.labels=c("ARMA","ARCH(8)","GARCH(1,1)","GARCH(1,2)","GARCH(2,1)","GARCH(2,2)","TGARCH(1,1)","TGARCH(1,2)","TGARCH(2,1)","TGARCH(2,2)"))

Estimation Results for Log Return of MSFT

	coredata	coredata	coredata	coredata	coredata	coredata	coredata	coredata	coredata	coredata
	ARIMA	GARCH	GARCH	GARCH	GARCH	GARCH	GARCH	GARCH	GARCH	GARCH
	ARMA	ARCH(8)	GARCH(1,1)	GARCH(1,2)	GARCH(2,1)	GARCH(2,2)	TGARCH(1,1)	TGARCH(1,2)	TGARCH(2,1)	TGARCH(2,2)
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
mu		0.0003***	0.0002***	0.0002***	0.0002***	0.0002***	0.0002***	0.0002***	0.0002***	0.0002***
		(0.0001)	(0.0001)	(0.0001)	(0.00005)	(0.00005)	(0.0001)	(0.0001)	(0.0001)	(0.0001)
ar1	-0.393**	0.286	0.317	0.312***	0.308***	0.286***	0.321	0.321	0.290***	0.321
	(0.183)			(0.004)	(0.006)	(0.005)			(0.005)
ar2	0.0001	0.991	0.869	0.878***	0.880***	0.920***	0.871	0.872	0.921***	0.872
	(0.160)			(0.011)	(0.004)	(0.006)			(0.005)
ar3	0.101	0.242***	0.460***	0.457***	0.457***	0.440***	0.456***	0.457***	0.436***	0.457***
	(0.168)	(0.062)	(0.014)	(0.015)	(0.004)	(0.005)	(0.014)	(0.014)	(0.005)	(0.014)
ar4	-0.601***	-0.794***	-0.898***	-0.896***	-0.897***	-0.895***	-0.903***	-0.902***	-0.901***	-0.902***
	(0.175)	(0.054)	(0.011)	(0.011)	(0.006)	(0.006)	(0.009)	(0.011)	(0.006)	(0.011)
ma1	0.355*	-0.312	-0.338***	-0.334***	-0.329***	-0.308***	-0.341***	-0.342***	-0.310***	-0.342***
	(0.184)		(0.007)	(0.008)	(0.001)	(0.003)	(0.008)	(0.009)	(0.003)	(0.009)
ma2	-0.062	-1.000	-0.891	-0.899***	-0.901***	-0.941***	-0.893	-0.893	-0.941***	-0.893
	(0.158)			(0.010)	(0.001)	(0.003)			(0.003)
ma3	-0.154	-0.254***	-0.455***	-0.452***	-0.454***	-0.439***	-0.452***	-0.453***	-0.436***	-0.453***
	(0.173)	(0.062)	(0.009)	(0.010)	(0.0002)	(0.001)	(0.008)	(0.007)	(0.001)	(0.007)
ma4	0.580***	0.813***	0.920***	0.920***	0.920***	0.918***	0.925***	0.924***	0.923***	0.924***
	(0.179)	(0.059)	(0.004)	(0.003)	(0.0002)	(0.001)	(0.002)	(0.005)	(0.001)	(0.005)
ma5	-0.028	0.008	0.002	0.002	0.003***	0.005**	0.003	0.003	0.006***	0.003
	(0.018)	(0.017)	(0.006)	(0.008)	(0.0004)	(0.002)	(0.005)	(0.004)	(0.002)	(0.004)
intercept	0.001***
	(0.0002)
omega		0.0002***	0.00000***	0.00001***	0.00001***	0.00001	0.00001***	0.00001***	0.00001***	0.00001***
		(0.00001)	(0.00000)	(0.00000)	(0.00000)		(0.00000)	(0.00000)	(0.00000)	(0.00000)
alpha1		0.160***	0.060***	0.075***	0.061***	0.074	0.060***	0.075***	0.061***	0.075***
		(0.015)	(0.006)	(0.008)	(0.012)		(0.006)	(0.008)	(0.006)	(0.008)
alpha2		0.124***			0.000	0.000			0.000	0.000
		(0.016)			(0.014)				(0.00001)	(0.00001)
alpha3		0.064***
		(0.012)
alpha4		0.115***
		(0.014)
alpha5		0.056***
		(0.012)
alpha6		0.030**
		(0.012)
alpha7		0.090***
		(0.014)
alpha8		0.054***
		(0.012)
gamma1							0.106***	0.108***	0.104***	0.108***
							(0.031)	(0.031)	(0.031)	(0.031)
gamma2									0.063	0.081
										(0.099)
beta1			0.932***	0.622***	0.931***	0.640	0.931***	0.635***	0.930***	0.635***
			(0.006)	(0.110)	(0.008)		(0.006)	(0.109)	(0.006)	(0.109)
beta2				0.292***		0.276		0.279***		0.279***
				(0.105)				(0.104)		(0.104)
Observations	7,584	7,584	7,584	7,584	7,584	7,584	7,584	7,584	7,584	7,584
Log Likelihood	17,810.100	-18,649.880	-18,834.410	-18,837.160	-18,836.210	-18,838.090	-18,840.540	-18,844.320	-18,841.140	-18,844.320
sigma2	0.001
Akaike Inf. Crit.	-35,598.200	-4.913	-4.963	-4.964	-4.964	-4.964	-4.965	-4.966	-4.964	-4.965
Bayesian Inf. Crit.		-4.896	-4.952	-4.951	-4.951	-4.950	-4.952	-4.952	-4.950	-4.949
skewness	-0.524	-0.185	-0.285	-0.284	-0.284	-0.298	-0.281	-0.274	-0.298	-0.274
kurtosis	14.473	6.424	6.759	6.705	6.773	6.779	6.606	6.513	6.721	6.513
Jarque-Bera	66,585.260 (p=0.000)	13,095.760 (p=0.000)	14,548.940 (p=0.000)	14,318.130 (p=0.000)	14,607.450 (p=0.000)	14,644.350 (p=0.000)	13,902.000 (p=0.000)	13,511.930 (p=0.000)	14,395.480 (p=0.000)	13,511.960 (p=0.000)
Q(10) resid	3.791 (p=0.957)	5.971 (p=0.818)	2.421 (p=0.992)	2.278 (p=0.994)	2.264 (p=0.994)	2.079 (p=0.996)	2.361 (p=0.993)	2.327 (p=0.994)	2.165 (p=0.995)	2.327 (p=0.994)
Q(20) resid	14.994 (p=0.777)	11.160 (p=0.942)	8.235 (p=0.991)	7.666 (p=0.994)	7.660 (p=0.994)	7.303 (p=0.996)	7.852 (p=0.993)	7.800 (p=0.994)	7.102 (p=0.997)	7.800 (p=0.994)
Q(10) resid sq.	740.784 (p=0.000)	6.798 (p=0.745)	1.919 (p=0.997)	1.693 (p=0.999)	1.765 (p=0.998)	1.778 (p=0.998)	2.059 (p=0.996)	1.905 (p=0.998)	1.970 (p=0.997)	1.905 (p=0.998)
Q(20) resid sq.	906.393 (p=0.000)	53.630 (p=0.0001)	8.380 (p=0.990)	7.862 (p=0.993)	8.410 (p=0.989)	7.953 (p=0.993)	8.643 (p=0.987)	8.143 (p=0.991)	8.593 (p=0.988)	8.143 (p=0.991)
sum phii	-0.894	0.725	0.748	0.751	0.748	0.752	0.746	0.748	0.746	0.748
sum (alphai+betai)		0.693	0.992	0.989	0.992	0.989	0.991	0.989	0.991	0.989
Note:	p<0.1; p<0.05; p<0.01

From our comparison, the best model based on AIC is actually TGARCH(1,2), we will this model as the base of our forecast.

(e) Forecasting

m11 <- ugarchfit(data=coredata(yt), spec=ugarchspec(mean.model=list(armaOrder=c(4,5), archm=TRUE, archpow=2, include.mean=TRUE), variance.model=list(model="sGARCH", garchOrder=c(1,2))), out.sample=200)
m11.fcst <- ugarchforecast(m11, n.roll=100, n.ahead=10)
plot(m11.fcst, which="all")