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")
