#Get data for the Microsoft stock price GOOG/NASDAQ_MSFT and construct the log return yt = ∆logMSFTt.
library("Quandl")
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library("tseries")
Quandl.auth(" 3HQPzReHC1WfzuSG41gy")
## Warning: 'Quandl.auth' is deprecated.
## Use 'Quandl.api_key' instead.
## See help("Deprecated")
msp<-Quandl("GOOG/NASDAQ_MSFT", type="zoo")
mspt<-msp$Close
lmspt<-log(mspt)
dlmspt <- diff(lmspt,1)
par(mfrow=c(1,2))
plot(lmspt,xlab="Time", ylab="lmspt")
plot(dlmspt,xlab="Time", ylab="dlmspt")

#a) Estimate an ARMA model: use ACF, PACF, Ljung-Box statistics, to #identify and then estimate an appropriate model for log return
plot(acf(coredata(dlmspt), type='correlation', lag=24, plot=FALSE), ylab="", main=expression(paste("ACF for dlmspt")))

acf(coredata(dlmspt),type="partial",lag=24,ylab="",main="PACF-dlmspt")

# estimate model -

m1 <- arima(dlmspt,order=c(2,0,2))
m1
## 
## Call:
## arima(x = dlmspt, order = c(2, 0, 2))
## 
## Coefficients:
##           ar1      ar2     ma1     ma2  intercept
##       -0.5344  -0.2307  0.4640  0.1463      9e-04
## s.e.   0.1690   0.1444  0.1717  0.1494      2e-04
## 
## sigma^2 estimated as 0.0005335:  log likelihood = 17803.19,  aic = -35594.37
#Check model m1 for adequacy
tsdiag(m1,gof.lag=24)

BIC(m1)
## [1] -35552.77
AIC(m1)
## [1] -35594.37
# try another model

m2 <- arima(dlmspt,order=c(4,0,4))
## Warning in arima(dlmspt, order = c(4, 0, 4)): possible convergence problem:
## optim gave code = 1
m2
## 
## Call:
## arima(x = dlmspt, order = c(4, 0, 4))
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ma1     ma2     ma3      ma4
##       0.6334  -0.4169  -0.0674  0.2443  -0.7075  0.4161  0.0978  -0.3218
## s.e.  0.2263   0.2177   0.2020  0.1604   0.2246  0.2281  0.2080   0.1559
##       intercept
##           8e-04
## s.e.      2e-04
## 
## sigma^2 estimated as 0.0005307:  log likelihood = 17815.3,  aic = -35610.6
#Check model m1 for adequacy
tsdiag(m2,gof.lag=24)

BIC(m2)
## [1] -35541.27
AIC(m2)
## [1] -35610.6

m2 is better no correlated resdulas and the value of BIC is lower than BIC for m1.

#(b) Use ACF, PACF, Ljung-Box statistics for squared residuals from (a) to identify and then estimate an ARCH(m) model with normal innovations
#combines the ARCH model for volatility with ARMA model for mean
par(mfrow=c(1,2))
acf(m2$residuals^2,type='correlation',na.action=na.pass,lag=24,main="ACF of Squared Residuals")
acf(m2$residuals^2,type='partial',na.action=na.pass,lag=24,main="PACF of Squared Residuals")

mm2 <- arima(m2$residuals^2, order=c(4,0,4))
tsdiag(m2,gof.lag=24)

library(fGarch)
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## 
## Rmetrics Package fBasics
## Analysing Markets and calculating Basic Statistics
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.org
n1 <- garchFit(~ arma(4,4) + garch(4,0), data = dlmspt, trace = FALSE) 
## Warning in arima(.series$x, order = c(u, 0, v), include.mean =
## include.mean): possible convergence problem: optim gave code = 1
## Warning in sqrt(diag(fit$cvar)): NaNs produced
n1
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 4) + garch(4, 0), data = dlmspt, 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 4) + garch(4, 0)
## <environment: 0x7fbf44e370d8>
##  [data = dlmspt]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  0.00021756   0.29746562   0.99999999   0.22258825  -0.78333927  
##         ma1          ma2          ma3          ma4        omega  
## -0.33068610  -0.99999999  -0.23359816   0.81271890   0.00023915  
##      alpha1       alpha2       alpha3       alpha4  
##  0.18626093   0.15802546   0.10569490   0.13068996  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      2.176e-04   6.038e-05    3.603 0.000315 ***
## ar1     2.975e-01          NA       NA       NA    
## ar2     1.000e+00          NA       NA       NA    
## ar3     2.226e-01   3.693e-02    6.027 1.67e-09 ***
## ar4    -7.833e-01   3.282e-03 -238.660  < 2e-16 ***
## ma1    -3.307e-01          NA       NA       NA    
## ma2    -1.000e+00          NA       NA       NA    
## ma3    -2.336e-01   3.776e-02   -6.187 6.12e-10 ***
## ma4     8.127e-01   6.830e-03  118.990  < 2e-16 ***
## omega   2.392e-04   7.457e-06   32.072  < 2e-16 ***
## alpha1  1.863e-01   1.625e-02   11.460  < 2e-16 ***
## alpha2  1.580e-01   1.760e-02    8.979  < 2e-16 ***
## alpha3  1.057e-01   1.344e-02    7.866 3.55e-15 ***
## alpha4  1.307e-01   1.463e-02    8.932  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18501.67    normalized:  2.439888 
## 
## Description:
##  Wed Mar 23 15:36:07 2016 by user:
#ARCH models - Checking Adequacy
par(mfrow=c(2,2), cex=0.6, mar=c(2,2,3,1))
plot(n1@residuals, type="l", xlab="", ylab="", main="residuals")
plot(n1@sigma.t, type="l", xlab="", ylab="", main="conditional standard deviation") 
plot(n1@residuals/n1@sigma.t, type="l", xlab="", ylab="", main="standardized residuals") 
plot((n1@residuals/n1@sigma.t)^2, type="l", xlab="", ylab="", main="squared standardized residuals")

#(c) Perform model checking for the ARCH model from (b): for both standardized residuals and squared standardized residuals plot ACF, examine Q(10) and Q(20) statistics. Also plot and comment on QQ plot and the results of Jarque-Bera test.
nu <- n1@residuals/n1@sigma.t
nu2 <- (n1@residuals/n1@sigma.t)^2
par(mfcol=c(2,2), cex=0.6, mar=c(2,2,3,1))
acf(nu, type="correlation", lag.max=24, ylab="", main="ACF for standardized residuals")
acf(nu, type="partial", lag.max=24, ylab="", main="PACF for standardized residuals")
acf(nu2, type="correlation", lag.max=24, ylab="", main="ACF for squared standardized residuals") 
acf(nu2, type="partial", lag.max=24, ylab="", main="PACF for squared standardized residuals")

Box.test( (n1@residuals/n1@sigma.t)^2, lag = 24, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  (n1@residuals/n1@sigma.t)^2
## X-squared = 158.22, df = 24, p-value < 2.2e-16
summary(n1)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(4, 4) + garch(4, 0), data = dlmspt, 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 4) + garch(4, 0)
## <environment: 0x7fbf44e370d8>
##  [data = dlmspt]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4  
##  0.00021756   0.29746562   0.99999999   0.22258825  -0.78333927  
##         ma1          ma2          ma3          ma4        omega  
## -0.33068610  -0.99999999  -0.23359816   0.81271890   0.00023915  
##      alpha1       alpha2       alpha3       alpha4  
##  0.18626093   0.15802546   0.10569490   0.13068996  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      2.176e-04   6.038e-05    3.603 0.000315 ***
## ar1     2.975e-01          NA       NA       NA    
## ar2     1.000e+00          NA       NA       NA    
## ar3     2.226e-01   3.693e-02    6.027 1.67e-09 ***
## ar4    -7.833e-01   3.282e-03 -238.660  < 2e-16 ***
## ma1    -3.307e-01          NA       NA       NA    
## ma2    -1.000e+00          NA       NA       NA    
## ma3    -2.336e-01   3.776e-02   -6.187 6.12e-10 ***
## ma4     8.127e-01   6.830e-03  118.990  < 2e-16 ***
## omega   2.392e-04   7.457e-06   32.072  < 2e-16 ***
## alpha1  1.863e-01   1.625e-02   11.460  < 2e-16 ***
## alpha2  1.580e-01   1.760e-02    8.979  < 2e-16 ***
## alpha3  1.057e-01   1.344e-02    7.866 3.55e-15 ***
## alpha4  1.307e-01   1.463e-02    8.932  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18501.67    normalized:  2.439888 
## 
## Description:
##  Wed Mar 23 15:36:07 2016 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  10642.81  0           
##  Shapiro-Wilk Test  R    W      NA        NA          
##  Ljung-Box Test     R    Q(10)  9.781869  0.4598348   
##  Ljung-Box Test     R    Q(15)  12.53819  0.6379196   
##  Ljung-Box Test     R    Q(20)  16.1841   0.7051396   
##  Ljung-Box Test     R^2  Q(10)  43.11444  4.744116e-06
##  Ljung-Box Test     R^2  Q(15)  72.67297  1.48602e-09 
##  Ljung-Box Test     R^2  Q(20)  139.3256  0           
##  LM Arch Test       R    TR^2   53.16328  3.852777e-07
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -4.876083 -4.863282 -4.876090 -4.871690
#d) Estimate a GARCH(m, s) model with Student-t innovations. Perform model checking.
garch22<-garchFit(~garch(2,2), data=coredata( dlmspt), trace=FALSE, cond.dist="sstd")
## Warning in sqrt(diag(fit$cvar)): NaNs produced
par(mfrow=c(2,2))
plot(garch22@residuals,type="l", xlab="",ylab="",main="Residuals")
plot(garch22@sigma.t,type="l", xlab="",ylab="",main="Conditional Standard Deviation")
plot(garch22@residuals/garch22@sigma.t, type="l", xlab="",ylab="",main="Standardized Residuals")
plot((garch22@residuals/garch22@sigma.t)^2, type="l", xlab="",ylab="",main="Squared Standardized Residuals")

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

#(e) Estimate a TGARCH(m, s) model, again perform model checking. What does this model tell us about
#the effects of negative vs positive shocks?#
tgarch<- garchFit(~garch(4,4), data=coredata(dlmspt), trace=FALSE, leverage=TRUE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
par(mfrow=c(2,2))
plot(tgarch@residuals,type="l", xlab="",ylab="",main="Residuals")
plot(tgarch@sigma.t,type="l", xlab="",ylab="",main="Conditional Standard Deviation")
plot(tgarch@residuals/tgarch@sigma.t, type="l", xlab="",ylab="",main="Standardized Residuals")
plot((tgarch@residuals/tgarch@sigma.t)^2, type="l", xlab="",ylab="",main="Squared Standardized Residuals")

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

#(f) Create a summary table using stargazerX showing the models side by side.
library("stargazer")
## 
## Please cite as:
##  Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2. http://CRAN.R-project.org/package=stargazer
stargazer(m1,n1,garch22,tgarch, align=TRUE, type="html", column.sep.width="0pt", report="vc*",
single.row=TRUE, header=FALSE, model.numbers=FALSE, dep.var.caption="", dep.var.labels="Microsoft Daily Returns", column.labels=c("ARMA(4,4)", "ARCH(4)", "GARCH(2,2)", "TGARCH(2,2)"))
## 
## <table style="text-align:center"><tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="2">Microsoft Daily Returns</td><td>coredata</td><td>coredata</td></tr>
## <tr><td style="text-align:left"></td><td><em>ARIMA</em></td><td><em>GARCH</em></td><td><em>GARCH</em></td><td><em>GARCH</em></td></tr>
## <tr><td style="text-align:left"></td><td>ARMA(4,4)</td><td>ARCH(4)</td><td>GARCH(2,2)</td><td>TGARCH(2,2)</td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">mu</td><td></td><td>0.0002<sup>***</sup></td><td>0.001<sup>***</sup></td><td>0.001<sup>***</sup></td></tr>
## <tr><td style="text-align:left">ar1</td><td>-0.534<sup>***</sup></td><td>0.297</td><td></td><td></td></tr>
## <tr><td style="text-align:left">ar2</td><td>-0.231</td><td>1.000</td><td></td><td></td></tr>
## <tr><td style="text-align:left">ar3</td><td></td><td>0.223<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left">ar4</td><td></td><td>-0.783<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left">ma1</td><td>0.464<sup>***</sup></td><td>-0.331</td><td></td><td></td></tr>
## <tr><td style="text-align:left">ma2</td><td>0.146</td><td>-1.000</td><td></td><td></td></tr>
## <tr><td style="text-align:left">intercept</td><td>0.001<sup>***</sup></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">ma3</td><td></td><td>-0.234<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left">ma4</td><td></td><td>0.813<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left">omega</td><td></td><td>0.0002<sup>***</sup></td><td>0.00000</td><td>0.00001</td></tr>
## <tr><td style="text-align:left">alpha1</td><td></td><td>0.186<sup>***</sup></td><td>0.077</td><td>0.077<sup>***</sup></td></tr>
## <tr><td style="text-align:left">alpha2</td><td></td><td>0.158<sup>***</sup></td><td>0.000</td><td>0.068</td></tr>
## <tr><td style="text-align:left">alpha3</td><td></td><td>0.106<sup>***</sup></td><td></td><td>0.000</td></tr>
## <tr><td style="text-align:left">alpha4</td><td></td><td>0.131<sup>***</sup></td><td></td><td>0.000</td></tr>
## <tr><td style="text-align:left">gamma1</td><td></td><td></td><td></td><td>0.219<sup>***</sup></td></tr>
## <tr><td style="text-align:left">gamma2</td><td></td><td></td><td></td><td>0.011</td></tr>
## <tr><td style="text-align:left">gamma3</td><td></td><td></td><td></td><td>0.140</td></tr>
## <tr><td style="text-align:left">gamma4</td><td></td><td></td><td></td><td>0.313</td></tr>
## <tr><td style="text-align:left">beta1</td><td></td><td></td><td>0.683</td><td>0.000</td></tr>
## <tr><td style="text-align:left">beta2</td><td></td><td></td><td>0.239</td><td>0.401</td></tr>
## <tr><td style="text-align:left">skew</td><td></td><td></td><td>1.041<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left">shape</td><td></td><td></td><td>4.833<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left">beta3</td><td></td><td></td><td></td><td>0.162</td></tr>
## <tr><td style="text-align:left">beta4</td><td></td><td></td><td></td><td>0.267<sup>***</sup></td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>7,583</td><td>7,583</td><td>7,583</td><td>7,583</td></tr>
## <tr><td style="text-align:left">Log Likelihood</td><td>17,803.190</td><td>-18,501.670</td><td>-19,266.870</td><td>-18,829.970</td></tr>
## <tr><td style="text-align:left">sigma<sup>2</sup></td><td>0.001</td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Akaike Inf. Crit.</td><td>-35,594.370</td><td>-4.876</td><td>-5.079</td><td>-4.963</td></tr>
## <tr><td style="text-align:left">Bayesian Inf. Crit.</td><td></td><td>-4.863</td><td>-5.072</td><td>-4.950</td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="4" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>