Homework2

(a) Build an AR model for the growth rate series. Perform model checking to validate the fitted model. Write down the model.

da1=diff((gdp.df[,2]))/gdp.df[-282,2]
pacf(da1,lag.max = 12)

Box.test(da1,lag=10,type = 'Ljung')

    Box-Ljung test

data:  da1
X-squared = 63.134, df = 10, p-value = 9.207e-10

t.test(da1)

    One Sample t-test

data:  da1
t = 13.772, df = 280, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.006698639 0.008932837
sample estimates:
  mean of x 
0.007815738 

We first build an AR(1) model roughly based on PACF

m0 <- arima(da1,order=c(1,0,0))
m0

Call:
arima(x = da1, order = c(1, 0, 0))

Coefficients:
         ar1  intercept
      0.3719     0.0078
s.e.  0.0553     0.0008

sigma^2 estimated as 7.765e-05:  log likelihood = 930.79,  aic = -1855.58

Model checking

tsdiag(m0,gof=20)

Refine the model

m1=ar(da1,method="mle") #Find the AR order
m1

Call:
ar(x = da1, method = "mle")

Coefficients:
      1        2        3  
 0.3455   0.1243  -0.0906  

Order selected 3  sigma^2 estimated as  7.633e-05

m1$aic

         0          1          2          3          4          5          6          7          8 
40.6706564  0.8178656  0.3617107  0.0000000  0.3438587  0.9452236  2.2069824  4.2006661  6.1562278 
         9         10         11         12 
 5.4566026  7.1833491  8.7169314  2.9566366 

m1$order

[1] 3

Model checking

Box.test(m1$resid,lag=10,type="Ljung")

    Box-Ljung test

data:  m1$resid
X-squared = 6.4222, df = 10, p-value = 0.7786

(pv1=1-pchisq(6.4412,10-3))

[1] 0.4892768

Fit an AR(3) model

(m2 <- arima(da1,order=c(3,0,0)))

Call:
arima(x = da1, order = c(3, 0, 0))

Coefficients:
         ar1     ar2      ar3  intercept
      0.3461  0.1255  -0.0930     0.0078
s.e.  0.0593  0.0624   0.0594     0.0008

sigma^2 estimated as 7.629e-05:  log likelihood = 933.26,  aic = -1856.51

Compute the constant term

phi = 1
for (i in 1:length(m2$coef)) {
if (i == length(m2$coef)) {
(phi = phi * m2$coef[i])
} else {
phi = phi-m2$coef[i]
}
}
print(phi)

        ar1 
0.004841478 

Model Checking

tsdiag(m2,gof=20)

Box.test(m2$resid,lag=10,type="Ljung")

    Box-Ljung test

data:  m2$resid
X-squared = 7.7724, df = 10, p-value = 0.6511

(pv2=1-pchisq(7.808,10-3))

[1] 0.3498285

Remove the insignificant coefficients. Since lag-3 coeff of m2 is insignificant at 5% level, we try AR(2)

(m3 <- arima(da1,order=c(2,0,0)))

Call:
arima(x = da1, order = c(2, 0, 0))

Coefficients:
         ar1     ar2  intercept
      0.3370  0.0940     0.0078
s.e.  0.0593  0.0593     0.0009

sigma^2 estimated as 7.696e-05:  log likelihood = 932.04,  aic = -1856.08

Model Checking

tsdiag(m3,gof=20)

Box.test(m3$resid,lag=10,type="Ljung")

    Box-Ljung test

data:  m3$resid
X-squared = 10.808, df = 10, p-value = 0.3726

(pv3 = 1-pchisq(10.746,10-2))

[1] 0.2165107

In the end, we find AR(3) model may be the best. The fitted model is \(r_{t}=0.0048+0.3461r_{t-1}+0.1255r_{t-2}-0.0930r_{t-3}+a_{t},\quad \{a_{t}\} \sim WN(0,7.629 \times 10^{-5})\)

(b) Does the model confirm the existence of business cycles? Why? (Hint: use the command polyroot to find roots of a polynomial.)

p1 <- c(1,-m2$coef[1:3])
(roots=polyroot(p1))

[1]  1.824805+1.15939i -2.300281+0.00000i  1.824805-1.15939i

Since the characteristic equation has complex solutions, the model confirms the existence of business cycles.

(c) Obtain 1-step to 8-step ahead point and 95% interval forecasts for the U.S. quarterly GDP growth rate at the forecast origin April 1, 2017 (the last data point).

(pre1=predict(m2,8))

$pred
Time Series:
Start = 282 
End = 289 
Frequency = 1 
[1] 0.007438792 0.008079120 0.007867656 0.007886415 0.007806813 0.007801287 0.007787639 0.007789626

$se
Time Series:
Start = 282 
End = 289 
Frequency = 1 
[1] 0.008734449 0.009242747 0.009487780 0.009492792 0.009493262 0.009494123 0.009494319 0.009494431

for (i in 1:length(pre1$pred)) {
  x1=pre1$pred[i]-1.96*pre1$se[i]
  x2=pre1$pred[i]+1.96*pre1$se[i]
  print(paste0(i,"-step ahead point and 95% interval forecast is","(",x1,", ",x2,")" ))
}

[1] "1-step ahead point and 95% interval forecast is(-0.00968072878869736, 0.0245583129248843)"
[1] "2-step ahead point and 95% interval forecast is(-0.0100366632540564, 0.0261949037065125)"
[1] "3-step ahead point and 95% interval forecast is(-0.0107283934291914, 0.026463704777231)"
[1] "4-step ahead point and 95% interval forecast is(-0.0107194567982441, 0.0264922873955126)"
[1] "5-step ahead point and 95% interval forecast is(-0.0107999795763257, 0.0264136064660958)"
[1] "6-step ahead point and 95% interval forecast is(-0.0108071940052144, 0.0264097671795465)"
[1] "7-step ahead point and 95% interval forecast is(-0.0108212266365856, 0.0263965048167583)"
[1] "8-step ahead point and 95% interval forecast is(-0.010819457971041, 0.0263987097547163)"

(d) Build an ARMA model for the log growth rate series. Perform model checking to validate the fitted model. Write down the model.

da2=diff(log(gdp.df[,2]))
acf(da2,lag.max = 12)

pacf(da2,lag.max = 12)

We first try ARMA(3,3)

m4 <- arima(da2,order=c(3,0,3))
m4

Call:
arima(x = da2, order = c(3, 0, 3))

Coefficients:
         ar1      ar2     ar3      ma1     ma2      ma3  intercept
      2.0804  -1.8559  0.5693  -1.7676  1.4043  -0.2943     0.0077
s.e.  0.1461   0.2191  0.1434   0.1646  0.2358   0.1564     0.0008

sigma^2 estimated as 7.167e-05:  log likelihood = 941.25,  aic = -1866.5

Model checking

tsdiag(m4,gof=20)

Box.test(m4$resid,lag=10,type="Ljung")

    Box-Ljung test

data:  m4$resid
X-squared = 5.223, df = 10, p-value = 0.8758

(pv4 = 1-pchisq(5.223,10-6))

[1] 0.2651719

Compute the constant term

ce=m4$coef
(phi0=ce[7]*(1-ce[1]-ce[2]-ce[3]))

  intercept 
0.001588681 

Remove the insignificant coefficients. Since coeff of ma3 is insignificant at 5% level, we try ARMA(3,2)

(m5 <- arima(da2,order=c(3,0,2)))

Call:
arima(x = da2, order = c(3, 0, 2))

Coefficients:
         ar1      ar2     ar3      ma1     ma2  intercept
      1.7143  -1.3237  0.2421  -1.3844  0.9182     0.0077
s.e.  0.1171   0.1346  0.0692   0.0897  0.0407     0.0007

sigma^2 estimated as 7.275e-05:  log likelihood = 939.62,  aic = -1865.25

Model checking

tsdiag(m5,gof=20)

Box.test(m5$resid,lag=10,type="Ljung")

    Box-Ljung test

data:  m5$resid
X-squared = 6.2192, df = 10, p-value = 0.7965

(pv5 = 1-pchisq(6.2192,10-5))

[1] 0.2854704

Since m4 has smaller AIC and larger p values for Ljung-box statistic, we choose ARMA(3,3). The fitted model is \(r_{t}=0.0016+2.080r_{t-1}-1.856r_{t-2}-0.569r_{t-3}+a_{t}+1.768a_{t-1}-1.404a_{t-2}+0.294a_{t-3},\quad \{a_{t}\} \sim WN(0,7.167 \times 10^{-5})\)

(e) Obtain 1-step to 8-step ahead point and 95% interval forecasts for the U.S. quarterly GDP growth rate at the forecast origin April 1, 2017 (the last data point).

(pre2=predict(m4,8))

$pred
Time Series:
Start = 282 
End = 289 
Frequency = 1 
[1] 0.005977463 0.006078627 0.006756573 0.007766862 0.008668049 0.009053858 0.008759187 0.007943227

$se
Time Series:
Start = 282 
End = 289 
Frequency = 1 
[1] 0.008466062 0.008870553 0.009029419 0.009076367 0.009081196 0.009082199 0.009087332 0.009089906

for (i in 1:length(pre1$pred)) {
  x1=pre1$pred[i]-1.96*pre1$se[i]
  x2=pre1$pred[i]+1.96*pre1$se[i]
  print(paste0(i,"-step ahead point and 95% interval forecast is","(",x1,", ",x2,")" ))
}

[1] "1-step ahead point and 95% interval forecast is(-0.00968072878869736, 0.0245583129248843)"
[1] "2-step ahead point and 95% interval forecast is(-0.0100366632540564, 0.0261949037065125)"
[1] "3-step ahead point and 95% interval forecast is(-0.0107283934291914, 0.026463704777231)"
[1] "4-step ahead point and 95% interval forecast is(-0.0107194567982441, 0.0264922873955126)"
[1] "5-step ahead point and 95% interval forecast is(-0.0107999795763257, 0.0264136064660958)"
[1] "6-step ahead point and 95% interval forecast is(-0.0108071940052144, 0.0264097671795465)"
[1] "7-step ahead point and 95% interval forecast is(-0.0108212266365856, 0.0263965048167583)"
[1] "8-step ahead point and 95% interval forecast is(-0.010819457971041, 0.0263987097547163)"

(a) Build a time series model for the log returns. Perform model checking to verify that the model is adequate, and write down the fitted model.

da4=log(da3+1)
pacf(da4) 

acf(da4) 

da5=diff(as.matrix(da4))
acf(da5)

da6=diff(da5,12)
acf(da6)

(m6 = arima(da4,order=c(1,0,1),seasonal=list(order=c(1,0,1),period=12)))

Call:
arima(x = da4, order = c(1, 0, 1), seasonal = list(order = c(1, 0, 1), period = 12))

Coefficients:

NaNs produced

         ar1     ma1    sar1     sma1  intercept
      0.0041  0.2143  0.9561  -0.8988     0.0104
s.e.     NaN     NaN     NaN      NaN     0.0052

sigma^2 estimated as 0.003605:  log likelihood = 1036.56,  aic = -2061.11

Model checking

tsdiag(m6,gof=24)

The fitted model is \((1-0.0041B)(1-0.9561^{12})r_{t}=(1-0.2143B)(1+0.8988B^{12})a_{t},\quad \{a_{t}\} \sim WN(0,0.003605)\)

(b) Create a January dummy variable by the command Jan = rep(c(1, rep(0, 11)), 62). Use the dummy variable as an explanatory variable and build a regression model with time series error for the log returns. Perform model checking to justify the model and write down the fitted model.

Jan = rep(c(1, rep(0, 11)), 62)
(m7 = arima(da4,order=c(0,0,1),xreg=Jan,include.mean=F))

Call:
arima(x = da4, order = c(0, 0, 1), xreg = Jan, include.mean = F)

Coefficients:
         ma1     Jan
      0.2372  0.0652
s.e.  0.0364  0.0073

sigma^2 estimated as 0.003479:  log likelihood = 1050.15,  aic = -2094.31

Model checking

tsdiag(m7,gof=24)

The model is \(r_{t}-0.0652Jan_{t}=a_{t}+0.2372a_{t-1},\quad \{a_{t}\} \sim WN(0,0.003479)\)

(c) Compare the two models in parts (a) and (b) using AIC criterion. Draw your conclusion about the seasonality of part (a).

Since AIC in (b) is smaller than (a), m7 is better than m6. And the log returns has a seasonality of 12 months.

(a) What model does the generated sequence rt1 simulate? What is its value at time 1 and 2?

rt1 is an ARMA(1,2) model Its value at time 1 and 2 is zero.

(b) Try to fit the series rt1 with a ARMA(1, 2) model and explain your finding.

(m8 <- arima(rt1,order=c(1,0,2)))

Call:
arima(x = rt1, order = c(1, 0, 2))

Coefficients:
         ar1      ma1     ma2  intercept
      0.5926  -0.0567  0.4448     1.7431
s.e.  0.0371   0.0377  0.0327     0.2209

sigma^2 estimated as 4.22:  log likelihood = -2139.41,  aic = 4288.81

The coeffient of ma2 is not equal to 2.

tsdiag(m8,gof=24)

(c) Replace the number of iteration 1000 to 100 and 10000, and repeat (b), (c)

at = rnorm(100)
rt1 = rep(0,100)
for (t in 3:100)
{
rt1[t]=0.5*rt1[t-1]+1+at[t]+2*at[t-2]
}
(m9 <- arima(rt1,order=c(1,0,2)))

Call:
arima(x = rt1, order = c(1, 0, 2))

Coefficients:
         ar1     ma1     ma2  intercept
      0.3809  0.1486  0.4978     0.9770
s.e.  0.1390  0.1185  0.1156     0.4497

sigma^2 estimated as 2.934:  log likelihood = -196.2,  aic = 402.4

The coeffient of ma2 is smaller.

at = rnorm(10000)
rt1 = rep(0,10000)
for (t in 3:10000)
{
rt1[t]=0.5*rt1[t-1]+1+at[t]+2*at[t-2]
}
(m10 <- arima(rt1,order=c(1,0,2)))

Call:
arima(x = rt1, order = c(1, 0, 2))

Coefficients:
         ar1      ma1     ma2  intercept
      0.5020  -0.0107  0.4904     1.9598
s.e.  0.0131   0.0124  0.0097     0.0600

sigma^2 estimated as 4.077:  log likelihood = -21216.53,  aic = 42443.06

The coeffient of ma2 approaches 2.

(d) Try to simulate another time series rt2 with the same model, using the same shock series at, but with different initial conditions rt2[1] = 20, rt2[2] = 10. What is the difference between rt1[1000] and rt2[1000]? Try to explain your result.

at = rnorm(1000)
rt1 = rep(0,1000)
for (t in 3:1000)
{
rt1[t]=0.5*rt1[t-1]+1+at[t]+2*at[t-2]
}
rt2 = rep(0,1000)
rt2[1] = 20
rt2[2] = 10
for (t in 3:1000)
{
rt2[t]=0.5*rt2[t-1]+1+at[t]+2*at[t-2]
}
rt1[1000]

[1] 5.834875

rt2[1000]

[1] 5.834875

rt1[1000] and rt2[1000] are the same. Because the first two time state have little effect.

(a) Compute the daily log returns of the S&P index. Is there any ARCH effect in the log returns?

[Hint: First fit an ARMA model for the returns, and then check the residuals of the fitted model for ARCH effect.]

da8=log(da7$sp+1)
acf(da8,lag.max = 12)

pacf(da8,lag.max = 12)

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

    Box-Ljung test

data:  da8
X-squared = 72.174, df = 12, p-value = 1.253e-10

t.test(da8)

    One Sample t-test

data:  da8
t = -0.45851, df = 2514, p-value = 0.6466
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.0006465277  0.0004014760
sample estimates:
    mean of x 
-0.0001225258 

(m10=arima(da8,c(2,0,2)))

Call:
arima(x = da8, order = c(2, 0, 2))

Coefficients:
          ar1      ar2     ma1     ma2  intercept
      -0.4727  -0.2164  0.3943  0.0965     -1e-04
s.e.   0.2111   0.2060  0.2166  0.2072      2e-04

sigma^2 estimated as 0.0001766:  log likelihood = 7297.98,  aic = -14583.96

require(aTSA)
arch.test(arima(m10$residuals,c(2,0,2)))

ARCH heteroscedasticity test for residuals 
alternative: heteroscedastic 

Portmanteau-Q test: 
     order   PQ p.value
[1,]     4  985       0
[2,]     8 2039       0
[3,]    12 3021       0
[4,]    16 3546       0
[5,]    20 4094       0
[6,]    24 4551       0
Lagrange-Multiplier test: 
     order   LM  p.value
[1,]     4 1367 0.00e+00
[2,]     8  444 0.00e+00
[3,]    12  270 0.00e+00
[4,]    16  188 0.00e+00
[5,]    20  146 0.00e+00
[6,]    24  113 7.16e-14

There is ARCH effect in the log returns.

(b) Fit a Gaussian ARMA-GARCH model for the log return series. Refine the model. Perform model checking, and write down the fitted model.

require(fGarch)
m11=garchFit(~arma(2,2)+garch(1,1),data=da8,trace=F)
summary(m11)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(2, 2) + garch(1, 1), data = da8, trace = F) 

Mean and Variance Equation:
 data ~ arma(2, 2) + garch(1, 1)
<environment: 0x0000022f8899d6d8>
 [data = da8]

Conditional Distribution:
 norm 

Coefficient(s):
         mu          ar1          ar2          ma1          ma2        omega       alpha1        beta1  
 9.4728e-05   5.6543e-01   9.6194e-02  -6.2678e-01  -1.0301e-01   1.0123e-06   7.1098e-02   9.2363e-01  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)    
mu      9.473e-05   6.431e-05    1.473 0.140734    
ar1     5.654e-01   4.382e-01    1.290 0.196920    
ar2     9.619e-02   3.381e-01    0.285 0.775989    
ma1    -6.268e-01   4.389e-01   -1.428 0.153297    
ma2    -1.030e-01   3.620e-01   -0.285 0.775949    
omega   1.012e-06   2.950e-07    3.431 0.000601 ***
alpha1  7.110e-02   9.107e-03    7.807 5.77e-15 ***
beta1   9.236e-01   9.721e-03   95.015  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 7863.915    normalized:  3.126805 

Description:
 Fri Jul 26 13:40:23 2019 by user: Pann 


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  258.0208  0           
 Shapiro-Wilk Test  R    W      0.9881416 1.296406e-13
 Ljung-Box Test     R    Q(10)  4.802282  0.9039883   
 Ljung-Box Test     R    Q(15)  14.57175  0.4826815   
 Ljung-Box Test     R    Q(20)  20.80817  0.4085      
 Ljung-Box Test     R^2  Q(10)  14.27224  0.1609322   
 Ljung-Box Test     R^2  Q(15)  17.97431  0.2640186   
 Ljung-Box Test     R^2  Q(20)  19.4759   0.4911131   
 LM Arch Test       R    TR^2   15.26362  0.2273357   

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-6.247249 -6.228704 -6.247269 -6.240518 

plot(m11)

Make a plot selection (or 0 to exit): 

 1:   Time Series                                 2:   Conditional SD                           
 3:   Series with 2 Conditional SD Superimposed   4:   ACF of Observations                      
 5:   ACF of Squared Observations                 6:   Cross Correlation                        
 7:   Residuals                                   8:   Conditional SDs                          
 9:   Standardized Residuals                     10:   ACF of Standardized Residuals            
11:   ACF of Squared Standardized Residuals      12:   Cross Correlation between r^2 and r      
13:   QQ-Plot of Standardized Residuals          

13

Make a plot selection (or 0 to exit): 

 1:   Time Series                                 2:   Conditional SD                           
 3:   Series with 2 Conditional SD Superimposed   4:   ACF of Observations                      
 5:   ACF of Squared Observations                 6:   Cross Correlation                        
 7:   Residuals                                   8:   Conditional SDs                          
 9:   Standardized Residuals                     10:   ACF of Standardized Residuals            
11:   ACF of Squared Standardized Residuals      12:   Cross Correlation between r^2 and r      
13:   QQ-Plot of Standardized Residuals          

0

Refine the model

require(fGarch)
m12=garchFit(~arma(1,1)+garch(2,1),data=da8,trace=F)
summary(m12)

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(1, 1) + garch(2, 1), data = da8, trace = F) 

Mean and Variance Equation:
 data ~ arma(1, 1) + garch(2, 1)
<environment: 0x0000022f86fd51d8>
 [data = da8]

Conditional Distribution:
 norm 

Coefficient(s):
         mu          ar1          ma1        omega       alpha1       alpha2        beta1  
 7.5797e-05   7.1977e-01  -7.8203e-01   1.3469e-06   3.0992e-03   8.2930e-02   9.0654e-01  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)    
mu      7.580e-05   4.444e-05    1.706 0.088087 .  
ar1     7.198e-01   9.122e-02    7.891 3.11e-15 ***
ma1    -7.820e-01   8.206e-02   -9.530  < 2e-16 ***
omega   1.347e-06   3.919e-07    3.437 0.000589 ***
alpha1  3.099e-03   1.656e-02    0.187 0.851501    
alpha2  8.293e-02   2.007e-02    4.132 3.59e-05 ***
beta1   9.065e-01   1.270e-02   71.365  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 7870.115    normalized:  3.129271 

Description:
 Fri Jul 26 13:40:31 2019 by user: Pann 


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  284.0684  0           
 Shapiro-Wilk Test  R    W      0.9872583 3.233526e-14
 Ljung-Box Test     R    Q(10)  4.524795  0.9205838   
 Ljung-Box Test     R    Q(15)  14.91451  0.4575922   
 Ljung-Box Test     R    Q(20)  21.09601  0.3915012   
 Ljung-Box Test     R^2  Q(10)  5.039164  0.8885465   
 Ljung-Box Test     R^2  Q(15)  8.726378  0.8913842   
 Ljung-Box Test     R^2  Q(20)  10.32324  0.9619249   
 LM Arch Test       R    TR^2   5.813987  0.9251668   

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-6.252975 -6.236748 -6.252990 -6.247085 

plot(m12)

Make a plot selection (or 0 to exit): 

 1:   Time Series                                 2:   Conditional SD                           
 3:   Series with 2 Conditional SD Superimposed   4:   ACF of Observations                      
 5:   ACF of Squared Observations                 6:   Cross Correlation                        
 7:   Residuals                                   8:   Conditional SDs                          
 9:   Standardized Residuals                     10:   ACF of Standardized Residuals            
11:   ACF of Squared Standardized Residuals      12:   Cross Correlation between r^2 and r      
13:   QQ-Plot of Standardized Residuals          

13

Make a plot selection (or 0 to exit): 

 1:   Time Series                                 2:   Conditional SD                           
 3:   Series with 2 Conditional SD Superimposed   4:   ACF of Observations                      
 5:   ACF of Squared Observations                 6:   Cross Correlation                        
 7:   Residuals                                   8:   Conditional SDs                          
 9:   Standardized Residuals                     10:   ACF of Standardized Residuals            
11:   ACF of Squared Standardized Residuals      12:   Cross Correlation between r^2 and r      
13:   QQ-Plot of Standardized Residuals          

0

Since m12 has smaller AIC, we choose m12. The fitted model is \(r_{t}=7.580\times10^{-5}+7.198\times10^{-1}r_{t-1}+a_{t}+7.820\times10^{-1}a_{t-1},\quad a_{t}=\sigma_{t}\epsilon_{t}, \quad \epsilon_{t} \sim N(0,1), \\ \sigma^{2}_{t}=1.347\times10^{-6}+3.099\times10^{-3}a^{2}_{t-1}+8.293\times10^{-2}a^{2}_{t-2}+9.065\times10^{-1}\sigma^{2}_{t-1}\)

(c) Compute 1-step to 4-step ahead forecasts of the log return and its volatility based on the fitted model.

predict(m12,4)