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
Quandl.api_key("4-KG5x_Vo7rXzmZNAHch")
library("tseries")
library("urca")
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
library(rugarch)
## Loading required package: parallel
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
ms <- Quandl("GOOG/NASDAQ_MSFT", type="zoo")
str(ms)
## 'zoo' series from 1986-03-13 to 2016-03-22
##   Data: num [1:7584, 1:5] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:5] "Open" "High" "Low" "Close" ...
##   Index:  Date[1:7584], format: "1986-03-13" "1986-03-14" "1986-03-17" "1986-03-18" ...
head(ms)
##            Open High Low Close   Volume
## 1986-03-13    0   NA  NA   0.1 30758400
## 1986-03-14    0   NA  NA   0.1  7852896
## 1986-03-17    0   NA  NA   0.1 33036768
## 1986-03-18    0   NA  NA   0.1 66403296
## 1986-03-19    0   NA  NA   0.1 47846304
## 1986-03-20    0   NA  NA   0.1 58377600
tail(ms)
##             Open  High   Low Close   Volume
## 2016-03-15 52.75 53.59 52.74 53.59 21104763
## 2016-03-16 53.45 54.60 53.40 54.35 31297892
## 2016-03-17 54.21 55.00 54.00 54.66 28186465
## 2016-03-18 54.92 54.97 53.44 53.49 67422879
## 2016-03-21 53.25 53.93 52.93 53.86 23807263
## 2016-03-22 53.61 54.25 53.46 54.07 23124143
ms2 <- rev(ms[,4])
lms<- log(ms2)
dlms <- diff(lms,1)
par(mfrow=c(1,2))
plot(lms,xlab="Periods", ylab="Log MS")
par(mfrow=c(1,2))

plot(ms, xlab ="Period", ylab= "closing price")

plot(dlms, xlab= "time", ylab= "%change")
acf(coredata(dlms), type = "correlation", lag = 24, main = "sample ACF")

acf(coredata(dlms), type = "partial", lag = 24, main = "sample PACF")
adf.test(dlms)
## Warning in adf.test(dlms): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlms
## Dickey-Fuller = -19.677, Lag order = 19, p-value = 0.01
## alternative hypothesis: stationary
m1 <- arima(coredata(dlms), order = c(4, 0, 0))
tsdiag(m1)

m1.LB <- Box.test(m1$residuals, lag=12, type="Ljung")
m1.LB
## 
##  Box-Ljung test
## 
## data:  m1$residuals
## X-squared = 8.7486, df = 12, p-value = 0.7242
m2 <- arima(coredata(dlms), order =c(5,0,0))
m2
## 
## Call:
## arima(x = coredata(dlms), order = c(5, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5  intercept
##       -0.0388  -0.0538  -0.0321  -0.0259  -0.0108     -8e-04
## s.e.   0.0115   0.0115   0.0115   0.0115   0.0115      2e-04
## 
## sigma^2 estimated as 0.0005352:  log likelihood = 17801.21,  aic = -35588.42
tsdiag(m2)

m2.LB <- Box.test(m2$residuals, lag=12, type="Ljung")
m2.LB
## 
##  Box-Ljung test
## 
## data:  m2$residuals
## X-squared = 7.9969, df = 12, p-value = 0.7854
acf(m2$residuals^2)
pacf(m2$residuals^2)

h1 <- garchFit(~ arma(3,2) + garch(8,0), data =dlms, trace = FALSE)
## Warning in arima(.series$x, order = c(u, 0, v), include.mean =
## include.mean): possible convergence problem: optim gave code = 1
h1
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(3, 2) + garch(8, 0), data = dlms, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(3, 2) + garch(8, 0)
## <environment: 0x7fd8f2011158>
##  [data = dlms]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ma1  
## -0.00061264  -0.15332450   0.50058179  -0.00238788   0.12323879  
##         ma2        omega       alpha1       alpha2       alpha3  
## -0.53260225   0.00015850   0.13999371   0.11919559   0.10665315  
##      alpha4       alpha5       alpha6       alpha7       alpha8  
##  0.10218279   0.07145601   0.04911015   0.08340365   0.06322144  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -6.126e-04   2.768e-04   -2.213 0.026869 *  
## ar1    -1.533e-01   2.853e-01   -0.537 0.591004    
## ar2     5.006e-01   1.931e-01    2.592 0.009550 ** 
## ar3    -2.388e-03   1.798e-02   -0.133 0.894351    
## ma1     1.232e-01   2.853e-01    0.432 0.665794    
## ma2    -5.326e-01   1.926e-01   -2.765 0.005684 ** 
## omega   1.585e-04   6.531e-06   24.269  < 2e-16 ***
## alpha1  1.400e-01   1.734e-02    8.075 6.66e-16 ***
## alpha2  1.192e-01   1.594e-02    7.478 7.53e-14 ***
## alpha3  1.067e-01   1.424e-02    7.489 6.93e-14 ***
## alpha4  1.022e-01   1.415e-02    7.223 5.08e-13 ***
## alpha5  7.146e-02   1.364e-02    5.239 1.61e-07 ***
## alpha6  4.911e-02   1.429e-02    3.437 0.000588 ***
## alpha7  8.340e-02   1.395e-02    5.978 2.26e-09 ***
## alpha8  6.322e-02   1.227e-02    5.152 2.58e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18717.35    normalized:  2.468331 
## 
## Description:
##  Wed Mar 23 17:03:07 2016 by user:
summary(h1)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(3, 2) + garch(8, 0), data = dlms, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(3, 2) + garch(8, 0)
## <environment: 0x7fd8f2011158>
##  [data = dlms]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ma1  
## -0.00061264  -0.15332450   0.50058179  -0.00238788   0.12323879  
##         ma2        omega       alpha1       alpha2       alpha3  
## -0.53260225   0.00015850   0.13999371   0.11919559   0.10665315  
##      alpha4       alpha5       alpha6       alpha7       alpha8  
##  0.10218279   0.07145601   0.04911015   0.08340365   0.06322144  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -6.126e-04   2.768e-04   -2.213 0.026869 *  
## ar1    -1.533e-01   2.853e-01   -0.537 0.591004    
## ar2     5.006e-01   1.931e-01    2.592 0.009550 ** 
## ar3    -2.388e-03   1.798e-02   -0.133 0.894351    
## ma1     1.232e-01   2.853e-01    0.432 0.665794    
## ma2    -5.326e-01   1.926e-01   -2.765 0.005684 ** 
## omega   1.585e-04   6.531e-06   24.269  < 2e-16 ***
## alpha1  1.400e-01   1.734e-02    8.075 6.66e-16 ***
## alpha2  1.192e-01   1.594e-02    7.478 7.53e-14 ***
## alpha3  1.067e-01   1.424e-02    7.489 6.93e-14 ***
## alpha4  1.022e-01   1.415e-02    7.223 5.08e-13 ***
## alpha5  7.146e-02   1.364e-02    5.239 1.61e-07 ***
## alpha6  4.911e-02   1.429e-02    3.437 0.000588 ***
## alpha7  8.340e-02   1.395e-02    5.978 2.26e-09 ***
## alpha8  6.322e-02   1.227e-02    5.152 2.58e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18717.35    normalized:  2.468331 
## 
## Description:
##  Wed Mar 23 17:03:07 2016 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  7421.354  0           
##  Shapiro-Wilk Test  R    W      NA        NA          
##  Ljung-Box Test     R    Q(10)  3.269518  0.9743512   
##  Ljung-Box Test     R    Q(15)  11.06242  0.7481573   
##  Ljung-Box Test     R    Q(20)  15.64134  0.7386107   
##  Ljung-Box Test     R^2  Q(10)  4.178894  0.9389161   
##  Ljung-Box Test     R^2  Q(15)  41.97729  0.0002265681
##  Ljung-Box Test     R^2  Q(20)  109.2729  2.664535e-14
##  LM Arch Test       R    TR^2   5.818726  0.9249426   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -4.932706 -4.918990 -4.932713 -4.927998
par(mfrow=c(2,2), cex = 0.5, mar = c(2,2,3,1))
plot(h1@residuals, xlab="", ylab="", main="Residuals")
plot(h1@sigma.t, xlab="", ylab="", main = "conditional standard deviation")
plot(h1@residuals/h1@sigma.t, xlab="", ylab="", main="Standarized residuals")
plot((h1@residuals/h1@sigma.t)^2, xlab="", ylab="", main = "squared standarized residuals")

acf1 <- h1@residuals/h1@sigma.t
acf2 <- (h1@residuals/h1@sigma.t)^2
par(mfcol=c(2,2), cex=0.5, mar=c(2,2,3,1))
acf(acf1, type="correlation", lag.max=18, ylab="", main="ACF for standarized residuals")
acf(acf1, type="partial", lag.max=18, ylab="", main="PACF for standarized residuals")
acf(acf2, type="correlation", lag.max = 18, ylab="", main="ACF for squared standarized residuals")
acf(acf2, type = "partial", lag.max=18, ylab="", main="PACF for squared standarized residuals")

sacf1 <- arima(acf1, order=c(10,0,0))
tsdiag(sacf1, gof.lag=18)

sacf1 <- arima(acf1, order=c(18,0,0))
tsdiag(sacf1,gof.lag=18)

sacf2 <- arima(acf2, order=c(10,0,0))
tsdiag(sacf2,gof.lag=18)

sacf2 <- arima(acf2, order=c(18,0,0))
tsdiag(sacf2,gof.lag=18)

arch1<- garchFit(~arma(3,2) + garch(8,0), data=dlms, trace=FALSE)
## Warning in arima(.series$x, order = c(u, 0, v), include.mean =
## include.mean): possible convergence problem: optim gave code = 1
par(mfrow=c(2,3), cex=0.5, mar=c(4,4,3,1))
plot(arch1, which=2)
plot(arch1, which=4)
plot(arch1, which=9)
plot(arch1, which=3)
plot(arch1, which=11)
plot(arch1, which=13)

skewness(dlms)
## [1] 0.4580773
## attr(,"method")
## [1] "moment"
kurtosis(dlms)
## [1] 14.45818
## attr(,"method")
## [1] "excess"
adf.test(dlms)
## Warning in adf.test(dlms): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlms
## Dickey-Fuller = -19.677, Lag order = 19, p-value = 0.01
## alternative hypothesis: stationary
normalTest(dlms, method="jb")
## 
## Title:
##  Jarque - Bera Normalality Test
## 
## Test Results:
##   STATISTIC:
##     X-squared: 66354.9472
##   P VALUE:
##     Asymptotic p Value: < 2.2e-16 
## 
## Description:
##  Wed Mar 23 17:04:46 2016 by user:
garch3 <- garchFit(~ garch(2,2), data =coredata(dlms), trace = FALSE, cond.dist="sstd")
## Warning in sqrt(diag(fit$cvar)): NaNs produced
garch3
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(2, 2), data = coredata(dlms), cond.dist = "sstd", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(2, 2)
## <environment: 0x7fd8f0b7ce28>
##  [data = coredata(dlms)]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu        omega       alpha1       alpha2        beta1  
## -8.1433e-04   3.7517e-06   9.1144e-02   1.0000e-08   7.0931e-01  
##       beta2         skew        shape  
##  1.9882e-01   9.4299e-01   5.0379e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -8.143e-04   1.899e-04   -4.287 1.81e-05 ***
## omega   3.752e-06          NA       NA       NA    
## alpha1  9.114e-02   1.014e-02    8.986  < 2e-16 ***
## alpha2  1.000e-08          NA       NA       NA    
## beta1   7.093e-01          NA       NA       NA    
## beta2   1.988e-01          NA       NA       NA    
## skew    9.430e-01   1.373e-02   68.679  < 2e-16 ***
## shape   5.038e+00   2.923e-01   17.233  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  19240.29    normalized:  2.537293 
## 
## Description:
##  Wed Mar 23 17:04:51 2016 by user:
par(mfrow=c(2,2))
tgarch<- garchFit(~garch(2,2), data=coredata(dlms), trace=FALSE, leverage=TRUE)
tgarch
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(2, 2), data = coredata(dlms), leverage = TRUE, 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(2, 2)
## <environment: 0x7fd8f141f990>
##  [data = coredata(dlms)]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1       alpha2       gamma1  
## -9.7839e-04   8.4850e-06   9.2631e-02   1.2973e-02   1.7406e-01  
##      gamma2        beta1        beta2  
## -1.0000e+00   5.5886e-01   3.1000e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -9.784e-04   2.005e-04   -4.881 1.06e-06 ***
## omega   8.485e-06   1.364e-06    6.220 4.96e-10 ***
## alpha1  9.263e-02   1.412e-02    6.562 5.31e-11 ***
## alpha2  1.297e-02   1.017e-02    1.276  0.20187    
## gamma1  1.741e-01   7.395e-02    2.354  0.01858 *  
## gamma2 -1.000e+00   7.266e-01   -1.376  0.16877    
## beta1   5.589e-01   1.159e-01    4.822 1.42e-06 ***
## beta2   3.100e-01   1.062e-01    2.919  0.00352 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18885.04    normalized:  2.490445 
## 
## Description:
##  Wed Mar 23 17:04:52 2016 by user:
m1
## 
## Call:
## arima(x = coredata(dlms), order = c(4, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3      ar4  intercept
##       -0.0385  -0.0534  -0.0315  -0.0255     -8e-04
## s.e.   0.0115   0.0115   0.0115   0.0115      2e-04
## 
## sigma^2 estimated as 0.0005352:  log likelihood = 17800.77,  aic = -35589.54
garch3
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(2, 2), data = coredata(dlms), cond.dist = "sstd", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(2, 2)
## <environment: 0x7fd8f0b7ce28>
##  [data = coredata(dlms)]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu        omega       alpha1       alpha2        beta1  
## -8.1433e-04   3.7517e-06   9.1144e-02   1.0000e-08   7.0931e-01  
##       beta2         skew        shape  
##  1.9882e-01   9.4299e-01   5.0379e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -8.143e-04   1.899e-04   -4.287 1.81e-05 ***
## omega   3.752e-06          NA       NA       NA    
## alpha1  9.114e-02   1.014e-02    8.986  < 2e-16 ***
## alpha2  1.000e-08          NA       NA       NA    
## beta1   7.093e-01          NA       NA       NA    
## beta2   1.988e-01          NA       NA       NA    
## skew    9.430e-01   1.373e-02   68.679  < 2e-16 ***
## shape   5.038e+00   2.923e-01   17.233  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  19240.29    normalized:  2.537293 
## 
## Description:
##  Wed Mar 23 17:04:51 2016 by user:
arch1
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(3, 2) + garch(8, 0), data = dlms, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(3, 2) + garch(8, 0)
## <environment: 0x7fd8f176e840>
##  [data = dlms]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ma1  
## -0.00061264  -0.15332450   0.50058179  -0.00238788   0.12323879  
##         ma2        omega       alpha1       alpha2       alpha3  
## -0.53260225   0.00015850   0.13999371   0.11919559   0.10665315  
##      alpha4       alpha5       alpha6       alpha7       alpha8  
##  0.10218279   0.07145601   0.04911015   0.08340365   0.06322144  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -6.126e-04   2.768e-04   -2.213 0.026869 *  
## ar1    -1.533e-01   2.853e-01   -0.537 0.591004    
## ar2     5.006e-01   1.931e-01    2.592 0.009550 ** 
## ar3    -2.388e-03   1.798e-02   -0.133 0.894351    
## ma1     1.232e-01   2.853e-01    0.432 0.665794    
## ma2    -5.326e-01   1.926e-01   -2.765 0.005684 ** 
## omega   1.585e-04   6.531e-06   24.269  < 2e-16 ***
## alpha1  1.400e-01   1.734e-02    8.075 6.66e-16 ***
## alpha2  1.192e-01   1.594e-02    7.478 7.53e-14 ***
## alpha3  1.067e-01   1.424e-02    7.489 6.93e-14 ***
## alpha4  1.022e-01   1.415e-02    7.223 5.08e-13 ***
## alpha5  7.146e-02   1.364e-02    5.239 1.61e-07 ***
## alpha6  4.911e-02   1.429e-02    3.437 0.000588 ***
## alpha7  8.340e-02   1.395e-02    5.978 2.26e-09 ***
## alpha8  6.322e-02   1.227e-02    5.152 2.58e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  18717.35    normalized:  2.468331 
## 
## Description:
##  Wed Mar 23 17:04:45 2016 by user:
tgarchspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), archm=TRUE, archpow.mean=TRUE), variance.model = list(model="sGARCH", garchOrder=c(2,2)))
## Warning: unidentified option(s) in mean.model:
##  archpow.mean
hf1<- ugarchfit(data=coredata(dlms), spec=tgarchspec,out.sample=200)
forecast <- ugarchforecast(hf1,n.roll = 100)
plot(forecast,which="all")