Importing microsoft stock price and calculate the log return:

mft <- Quandl("GOOG/NASDAQ_MSFT", type="zoo")
lmft<- diff(log(mft[,4]))

Plotting the series at level and log transformation:

par(mfrow=c(1,2), cex=0.8, mar=c(2,2,3,1))
plot(mft[,4], xlab="", ylab="Dollars", main=" Closing Price")
plot(lmft, xlab="", ylab="Log", main="Log Return of Microsoft")

Part a: ARMA

ACF and PACF:

From the ACF and PACF, we will use ARMA(4,0,2) specification:

arma <- arima(coredata(lmft),order=c(4,0,2))
arma
## 
## Call:
## arima(x = coredata(lmft), order = c(4, 0, 2))
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ma1      ma2  intercept
##       -0.0195  -0.0265  -0.0298  -0.0234  -0.0191  -0.0262      8e-04
## s.e.      NaN      NaN   0.0117   0.0117      NaN      NaN      2e-04
## 
## sigma^2 estimated as 0.0005351:  log likelihood = 17806.47,  aic = -35596.94
tsdiag(arma,gof.lag=12)

Based on the plots above, ARIMA(4,0,2)specification is adequate. Next we will test for Ljung-Box statistics at lag 8 (based on the number of the observation):

Box.test(coredata(arma$residuals), lag=8, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  coredata(arma$residuals)
## X-squared = 4.5946, df = 8, p-value = 0.7999

The result confirmed our previous adequacy check.

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

The Ljung-Box test of the residuals of ARMA model previously estimated shows that there is no serial dependence in the residuals, hence, we can proceed to estimating the joint ARMA-ARCH model. Next we will use the ACF and PACF of the residuals to identify the order of \(m\) for our ARCH(\(m\)) model.

From the PACF big spikes from lag 1 to lag 8, we have an ARCH(8) model. We estimate the ARCH(8) model as follow:

armagarch <- garchFit(~ arma(4,2) + garch(8,0), data = coredata(lmft), trace = FALSE)
armagarch
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 2) + garch(8, 0), data = coredata(lmft), 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 2) + garch(8, 0)
## <environment: 0x000000001f248950>
##  [data = coredata(lmft)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  0.00031526   0.97608725  -0.29658428  -0.00302696  -0.01524468  
##         ma1          ma2        omega       alpha1       alpha2  
## -0.99999999   0.29389402   0.00018142   0.16111606   0.12060745  
##      alpha3       alpha4       alpha5       alpha6       alpha7  
##  0.06216388   0.11268369   0.05480254   0.02868124   0.09437387  
##      alpha8  
##  0.05372201  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      3.153e-04   5.686e-05    5.544 2.95e-08 ***
## ar1     9.761e-01          NA       NA       NA    
## ar2    -2.966e-01          NA       NA       NA    
## ar3    -3.027e-03   1.753e-02   -0.173   0.8629    
## ar4    -1.524e-02   1.443e-02   -1.057   0.2907    
## ma1    -1.000e+00          NA       NA       NA    
## ma2     2.939e-01          NA       NA       NA    
## omega   1.814e-04   7.054e-06   25.718  < 2e-16 ***
## alpha1  1.611e-01   1.501e-02   10.736  < 2e-16 ***
## alpha2  1.206e-01   1.607e-02    7.504 6.17e-14 ***
## alpha3  6.216e-02   1.212e-02    5.130 2.90e-07 ***
## alpha4  1.127e-01   1.439e-02    7.829 4.88e-15 ***
## alpha5  5.480e-02   1.145e-02    4.788 1.68e-06 ***
## alpha6  2.868e-02   1.163e-02    2.466   0.0136 *  
## alpha7  9.437e-02   1.446e-02    6.526 6.75e-11 ***
## alpha8  5.372e-02   1.175e-02    4.572 4.82e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18640.93    normalized:  2.457604 
## 
## Description:
##  Fri Mar 25 23:11:25 2016 by user: evanu

Part c: Checking for model adequacy

Constructing the standardized residuals from previous ARCH model:

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

resid.std <- armagarch@residuals/armagarch@sigma.t
sqr.resid.std <- (armagarch@residuals/armagarch@sigma.t)^2

plot(resid.std, type="l", xlab="", ylab="", main="standardized residuals")
plot(sqr.resid.std, type="l", xlab="", ylab="", main="squared standardized residuals")

Plot the ACF and PACF:

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

Ljung Box test for Q(10) and Q(20)

Box.test(sqr.resid.std, lag = 10, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sqr.resid.std
## X-squared = 7.1123, df = 10, p-value = 0.7148
Box.test(sqr.resid.std, lag = 20, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sqr.resid.std
## X-squared = 52.513, df = 20, p-value = 9.582e-05

From the test above, we can see that our model is adequate, indicated by no correlation at Q(10) eventhough the result is borderline for Q(20).

QQplot:

plot(armagarch, which=13)

Part d: GARCH Model

Estimating our GARCH(1,1) model:

garch11 <- garchFit(~ arma(4,2) + garch(1,1), data = coredata(lmft), trace = FALSE, cond.dist="sstd")
summary(garch11)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 2) + garch(1, 1), data = coredata(lmft), 
##     cond.dist = "sstd", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 2) + garch(1, 1)
## <environment: 0x000000001a3e0e68>
##  [data = coredata(lmft)]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  2.7434e-04   9.6738e-01  -3.1958e-01   1.0624e-02  -1.8846e-02  
##         ma1          ma2        omega       alpha1        beta1  
## -1.0000e+00   3.2941e-01   2.5760e-06   6.5113e-02   9.3392e-01  
##        skew        shape  
##  1.0381e+00   4.8167e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      2.743e-04   8.848e-05    3.101 0.001932 ** 
## ar1     9.674e-01          NA       NA       NA    
## ar2    -3.196e-01          NA       NA       NA    
## ar3     1.062e-02   1.545e-02    0.688 0.491647    
## ar4    -1.885e-02   8.341e-03   -2.260 0.023852 *  
## ma1    -1.000e+00          NA       NA       NA    
## ma2     3.294e-01          NA       NA       NA    
## omega   2.576e-06   7.219e-07    3.568 0.000359 ***
## alpha1  6.511e-02   8.642e-03    7.534 4.91e-14 ***
## beta1   9.339e-01   8.311e-03  112.368  < 2e-16 ***
## skew    1.038e+00   1.490e-02   69.678  < 2e-16 ***
## shape   4.817e+00   2.680e-01   17.974  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  19279.53    normalized:  2.541797 
## 
## Description:
##  Fri Mar 25 23:12:25 2016 by user: evanu 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  23767.67  0        
##  Shapiro-Wilk Test  R    W      NA        NA       
##  Ljung-Box Test     R    Q(10)  4.950314  0.8944722
##  Ljung-Box Test     R    Q(15)  11.08524  0.7465284
##  Ljung-Box Test     R    Q(20)  14.15302  0.8226536
##  Ljung-Box Test     R^2  Q(10)  1.853821  0.9973427
##  Ljung-Box Test     R^2  Q(15)  7.080892  0.955366 
##  Ljung-Box Test     R^2  Q(20)  8.378297  0.9890542
##  LM Arch Test       R    TR^2   2.859808  0.9964554
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -5.080431 -5.069461 -5.080436 -5.076666

Constructing the standardized residuals from previous GARCH model:

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

QQplot:

plot(garch11, which=13)

Part e: TGARCH Model

Estimating our TGARCH(1,1) model:

tgarch1 <- garchFit(~ arma(4,2) + garch(1,1), data=coredata(lmft), trace=FALSE, leverage=TRUE)
summary(tgarch1)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 2) + garch(1, 1), data = coredata(lmft), 
##     leverage = TRUE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 2) + garch(1, 1)
## <environment: 0x000000001fac5e48>
##  [data = coredata(lmft)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  2.3004e-04   9.8082e-01  -3.1294e-01   8.8886e-03  -1.8088e-02  
##         ma1          ma2        omega       alpha1       gamma1  
## -1.0000e+00   3.0356e-01   5.1905e-06   5.9204e-02   9.4139e-02  
##       beta1  
##  9.3182e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      2.300e-04   6.796e-05    3.385 0.000712 ***
## ar1     9.808e-01          NA       NA       NA    
## ar2    -3.129e-01          NA       NA       NA    
## ar3     8.889e-03   1.528e-02    0.582 0.560711    
## ar4    -1.809e-02   1.229e-02   -1.472 0.140939    
## ma1    -1.000e+00          NA       NA       NA    
## ma2     3.036e-01          NA       NA       NA    
## omega   5.190e-06   7.729e-07    6.715 1.88e-11 ***
## alpha1  5.920e-02   5.501e-03   10.762  < 2e-16 ***
## gamma1  9.414e-02   3.057e-02    3.080 0.002072 ** 
## beta1   9.318e-01   6.126e-03  152.098  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18830.92    normalized:  2.482652 
## 
## Description:
##  Fri Mar 25 23:12:34 2016 by user: evanu 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value  
##  Jarque-Bera Test   R    Chi^2  14915.45  0        
##  Shapiro-Wilk Test  R    W      NA        NA       
##  Ljung-Box Test     R    Q(10)  3.025138  0.9808263
##  Ljung-Box Test     R    Q(15)  8.803168  0.8875821
##  Ljung-Box Test     R    Q(20)  11.39793  0.9352432
##  Ljung-Box Test     R^2  Q(10)  2.090336  0.9955997
##  Ljung-Box Test     R^2  Q(15)  7.163728  0.9529461
##  Ljung-Box Test     R^2  Q(20)  8.74592   0.9856676
##  LM Arch Test       R    TR^2   2.875479  0.9963611
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -4.962404 -4.952348 -4.962408 -4.958953

Constructing the standardized residuals from previous TGARCH model:

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

plot(tgarch1, which=13)

Part f: Stargazer

Loading the source:

setwd("C:/Users/evanu/OneDrive/Time Series/Homework 4")
source("stargazerX.r")
source("stargazerX-internal.r")

We compile our stargazer table using the following command:

stargazerX(arma,armagarch,garch11,tgarch1, type="html", title="", 
           align=TRUE, header=FALSE, model.numbers=TRUE,
           dep.var.caption="", 
           column.labels=c("ARMA","ARMA-ARCH(8)","GARCH(1,1)","TGARCH"))
coredata coredata coredata coredata
ARIMA GARCH GARCH GARCH
ARMA ARMA-ARCH(8) GARCH(1,1) TGARCH
(1) (2) (3) (4)
mu 0.0003*** 0.0003*** 0.0002***
(0.0001) (0.0001) (0.0001)
ar1 -0.020 0.976 0.967 0.981
ar2 -0.026 -0.297 -0.320 -0.313
ar3 -0.030** -0.003 0.011 0.009
(0.012) (0.018) (0.015) (0.015)
ar4 -0.023** -0.015 -0.019** -0.018
(0.012) (0.014) (0.008) (0.012)
ma1 -0.019 -1.000 -1.000 -1.000
ma2 -0.026 0.294 0.329 0.304
intercept 0.001***
(0.0002)
omega 0.0002*** 0.00000*** 0.00001***
(0.00001) (0.00000) (0.00000)
alpha1 0.161*** 0.065*** 0.059***
(0.015) (0.009) (0.006)
alpha2 0.121***
(0.016)
alpha3 0.062***
(0.012)
alpha4 0.113***
(0.014)
alpha5 0.055***
(0.011)
alpha6 0.029**
(0.012)
alpha7 0.094***
(0.014)
alpha8 0.054***
(0.012)
gamma1 0.094***
(0.031)
beta1 0.934*** 0.932***
(0.008) (0.006)
skew 1.038***
(0.015)
shape 4.817***
(0.268)
Observations 7,585 7,585 7,585 7,585
Log Likelihood 17,806.470 -18,640.930 -19,279.530 -18,830.920
sigma2 0.001
Akaike Inf. Crit. -35,596.940 -4.911 -5.080 -4.962
Bayesian Inf. Crit. -4.896 -5.069 -4.952
skewness -0.533 -0.211 -0.395 -0.312
kurtosis 14.661 6.637 8.633 6.839
Jarque-Bera 68,330.990 (p=0.000) 13,990.060 (p=0.000) 23,767.670 (p=0.000) 14,915.450 (p=0.000)
Q(10) resid 8.231 (p=0.607) 3.759 (p=0.958) 4.950 (p=0.895) 3.025 (p=0.981)
Q(20) resid 20.042 (p=0.456) 10.546 (p=0.958) 14.153 (p=0.823) 11.398 (p=0.936)
Q(10) resid sq. 768.211 (p=0.000) 7.112 (p=0.715) 1.854 (p=0.998) 2.090 (p=0.996)
Q(20) resid sq. 936.947 (p=0.000) 52.513 (p=0.0001) 8.378 (p=0.990) 8.746 (p=0.986)
sum phii -0.099 0.661 0.640 0.659
sum (alphai+betai) 0.688 0.999 0.991
Note: p<0.1; p<0.05; p<0.01

By comparing all of our estimated model, the best model is GARCH(1,1), hence, we will use this model in the next part to be the base of our forecast.

Part g: Forecasting using rugarch Packages

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