Untitled
King Harbor Wealth Management
4/25/2021
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(quantmod)## Warning: package 'quantmod' was built under R version 4.0.4## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Warning: package 'TTR' was built under R version 4.0.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries)## Warning: package 'tseries' was built under R version 4.0.4library(fGarch)## Warning: package 'fGarch' was built under R version 4.0.4## Loading required package: timeDate## Loading required package: timeSeries## Warning: package 'timeSeries' was built under R version 4.0.2##
## Attaching package: 'timeSeries'## The following object is masked from 'package:zoo':
##
## time<-## Loading required package: fBasics## Warning: package 'fBasics' was built under R version 4.0.2##
## Attaching package: 'fBasics'## The following object is masked from 'package:TTR':
##
## volatilitylibrary(mgcv)## Loading required package: nlme## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.library(rugarch)## Warning: package 'rugarch' was built under R version 4.0.4## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigmarm(list = ls())
getwd()## [1] "C:/Users/perei/Dropbox/Worksheets/Lawrence/Python/Time Series Analysis - ISYE-6402-OANO01/Final"setwd("C:/Users/perei/Dropbox/Worksheets/Lawrence/Python/Time Series Analysis - ISYE-6402-OANO01/Final")
#data=read.csv("EGDailyVolume.csv",header=T)
gnp = read.csv("GNP.csv",header=T)
gnp=gnp[,2]
gnpgr=diff(log(gnp))
ts.plot(gnpgr,ylab="the growth rate of GNP")
n=length(gnpgr)
gnpgr.test = gnpgr[200:n]
gnpgr.train = gnpgr[-c(200:n)]
#Q1
#ARIMA order
test_modelA <- function(p,d,q){
mod = arima(gnpgr.train, order=c(p,d,q), method="ML")
current.aic = AIC(mod)
df = data.frame(p,d,q,current.aic)
names(df) <- c("p","d","q","AIC")
print(paste(p,d,q,current.aic,sep=" "))
return(df)
}
orders = data.frame(Inf,Inf,Inf,Inf)
names(orders) <- c("p","d","q","AIC")
for (p in 0:4){
for (d in 0:1){
for (q in 0:4) {
possibleError <- tryCatch(
orders<-rbind(orders,test_modelA(p,d,q)),
error=function(e) e
)
if(inherits(possibleError, "error")) next
}
}
}## [1] "0 0 0 -1204.51843294805"
## [1] "0 0 1 -1234.84861443414"
## [1] "0 0 2 -1251.12112705159"
## [1] "0 0 3 -1252.0458228431"
## [1] "0 0 4 -1252.83440072015"
## [1] "0 1 0 -1184.63877912904"
## [1] "0 1 1 -1211.17756004753"
## [1] "0 1 2 -1228.23391009149"
## [1] "0 1 3 -1241.20864741795"
## [1] "0 1 4 -1241.16033089472"
## [1] "1 0 0 -1249.81662882813"
## [1] "1 0 1 -1249.10200102136"
## [1] "1 0 2 -1252.13454907672"
## [1] "1 0 3 -1250.6820014958"
## [1] "1 0 4 -1253.24480809898"
## [1] "1 1 0 -1208.10141386723"
## [1] "1 1 1 -1238.98074388342"
## [1] "1 1 2 -1237.93563721728"
## [1] "1 1 3 -1241.06467552833"
## [1] "1 1 4 -1239.57115528764"
## [1] "2 0 0 -1250.10663641891"
## [1] "2 0 1 -1251.23709976674"
## [1] "2 0 2 -1255.38430710539"
## [1] "2 0 3 -1254.49158939393"
## [1] "2 0 4 -1254.04057799107"
## [1] "2 1 0 -1206.4265099929"
## [1] "2 1 1 -1238.95385186291"
## [1] "2 1 2 -1240.04040906855"
## [1] "2 1 3 -1238.89465784869"
## [1] "2 1 4 -1237.37152317369"
## [1] "3 0 0 -1253.75266191653"
## [1] "3 0 1 -1253.36279051425"
## [1] "3 0 2 -1257.52438731093"## Warning in arima(gnpgr.train, order = c(p, d, q), method = "ML"): possible
## convergence problem: optim gave code = 1## [1] "3 0 3 -1252.59886816811"
## [1] "3 0 4 -1256.73025533805"
## [1] "3 1 0 -1208.29219032123"## Warning in log(s2): NaNs produced## Warning in log(s2): NaNs produced## [1] "3 1 1 -1243.76781773491"
## [1] "3 1 2 -1249.04161976605"
## [1] "3 1 3 -1247.43326978656"
## [1] "3 1 4 -1246.56233364531"
## [1] "4 0 0 -1252.87319979692"
## [1] "4 0 1 -1251.52765268104"
## [1] "4 0 2 -1249.48742591571"
## [1] "4 0 3 -1255.78313399774"
## [1] "4 0 4 -1254.0955416042"
## [1] "4 1 0 -1208.01078944559"
## [1] "4 1 1 -1244.43614651682"
## [1] "4 1 2 -1247.3691794681"
## [1] "4 1 3 -1246.00778167601"
## [1] "4 1 4 -1237.96438895161"orders <- orders[order(-orders$AIC),]
tail(orders)## p d q AIC
## 46 4 0 4 -1254.096
## 25 2 0 3 -1254.492
## 24 2 0 2 -1255.384
## 45 4 0 3 -1255.783
## 36 3 0 4 -1256.730
## 34 3 0 2 -1257.524#3,0,2
final.arima = arima(gnpgr.train, order=c(3,0,2))
## Residual Analysis
resids = resid(final.arima)
par(mfrow=c(1,2))
acf(resids,main="Residuals of ARIMA Fit")
acf(resids^2,main="Squared Residuals of ARIMA Fit")
# for serial correlation
Box.test(resids,lag=6,type='Ljung',fitdf=5)##
## Box-Ljung test
##
## data: resids
## X-squared = 7.1173, df = 1, p-value = 0.007634# for heteroscedasticity in residuals (ARCH effect)
Box.test((resids)^2,lag=6,type='Ljung',fitdf=5)##
## Box-Ljung test
##
## data: (resids)^2
## X-squared = 19.688, df = 1, p-value = 9.117e-06#Q2
#ARIMA-GARCH GARCH order
#GARCH update
test_modelAGG <- function(m,n){
spec = ugarchspec(variance.model=list(garchOrder=c(m,n)),
mean.model=list(armaOrder=c(3,2),
include.mean=T), distribution.model="std")
fit = ugarchfit(spec, gnpgr.train, solver = 'hybrid')
current.bic = infocriteria(fit)[2]
df = data.frame(m,n,current.bic)
names(df) <- c("m","n","BIC")
print(paste(m,n,current.bic,sep=" "))
return(df)
}
orders = data.frame(Inf,Inf,Inf)
names(orders) <- c("m","n","BIC")
for (m in 0:2){
for (n in 0:2){
possibleError <- tryCatch(
orders<-rbind(orders,test_modelAGG(m,n)),
error=function(e) e
)
if(inherits(possibleError, "error")) next
}
}## [1] "0 0 -6.25569047025327"
## [1] "0 1 -6.23668744372467"
## [1] "0 2 -6.20980771122167"
## [1] "1 0 -6.24830850208647"
## [1] "1 1 -6.30713036002182"
## [1] "1 2 -6.280530728305"
## [1] "2 0 -6.23019475967242"
## [1] "2 1 -6.28065355666198"
## [1] "2 2 -6.25581336783202"orders <- orders[order(-orders$BIC),]
tail(orders)## m n BIC
## 5 1 0 -6.248309
## 2 0 0 -6.255690
## 10 2 2 -6.255813
## 7 1 2 -6.280531
## 9 2 1 -6.280654
## 6 1 1 -6.307130#1,1
#ARMA update
#ARIMA-GARCH ARIMA order
test_modelAGA <- function(p,q){
spec = ugarchspec(variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(p,q),
include.mean=T), distribution.model="std")
fit = ugarchfit(spec, gnpgr.train, solver = 'hybrid')
current.bic = infocriteria(fit)[2]
df = data.frame(p,q,current.bic)
names(df) <- c("p","q","BIC")
print(paste(p,q,current.bic,sep=" "))
return(df)
}
orders = data.frame(Inf,Inf,Inf)
names(orders) <- c("p","q","BIC")
for (p in 0:4){
for (q in 0:4){
possibleError <- tryCatch(
orders<-rbind(orders,test_modelAGA(p,q)),
error=function(e) e
)
if(inherits(possibleError, "error")) next
}
}## [1] "0 0 -6.1672777247186"
## [1] "0 1 -6.2740663735071"
## [1] "0 2 -6.35366421650449"
## [1] "0 3 -6.3476638691336"
## [1] "0 4 -6.3261560665703"
## [1] "1 0 -6.35143598926223"
## [1] "1 1 -6.33941598638969"
## [1] "1 2 -6.34625205576785"
## [1] "1 3 -6.32398938808612"
## [1] "1 4 -6.30667458236625"
## [1] "2 0 -6.35253720301128"
## [1] "2 1 -6.34053171539081"
## [1] "2 2 -6.32854175278947"
## [1] "2 3 -6.30224259917422"
## [1] "2 4 -6.2790633591618"
## [1] "3 0 -6.35372103140011"
## [1] "3 1 -6.33020830162369"
## [1] "3 2 -6.30713036002182"
## [1] "3 3 -6.29612476488999"
## [1] "3 4 -6.31546113512338"
## [1] "4 0 -6.32462773535165"
## [1] "4 1 -6.298371518636"
## [1] "4 2 -6.28907916613839"
## [1] "4 3 -6.26916549562144"
## [1] "4 4 -6.23687600721414"orders <- orders[order(-orders$BIC),]
tail(orders)## p q BIC
## 9 1 2 -6.346252
## 5 0 3 -6.347664
## 7 1 0 -6.351436
## 12 2 0 -6.352537
## 4 0 2 -6.353664
## 17 3 0 -6.353721#3,0
#Final Garch Order
#GARCH update
test_modelAGG <- function(m,n){
spec = ugarchspec(variance.model=list(garchOrder=c(m,n)),
mean.model=list(armaOrder=c(3,0),
include.mean=T), distribution.model="std")
fit = ugarchfit(spec, gnpgr.train, solver = 'hybrid')
current.bic = infocriteria(fit)[2]
df = data.frame(m,n,current.bic)
names(df) <- c("m","n","BIC")
print(paste(m,n,current.bic,sep=" "))
return(df)
}
orders = data.frame(Inf,Inf,Inf)
names(orders) <- c("m","n","BIC")
for (m in 0:2){
for (n in 0:2){
possibleError <- tryCatch(
orders<-rbind(orders,test_modelAGG(m,n)),
error=function(e) e
)
if(inherits(possibleError, "error")) next
}
}## [1] "0 0 -6.27381336408152"
## [1] "0 1 -6.27665887319072"
## [1] "0 2 -6.24980398934888"
## [1] "1 0 -6.29525141242357"
## [1] "1 1 -6.35372103140011"
## [1] "1 2 -6.32717933734311"
## [1] "2 0 -6.30857412253723"
## [1] "2 1 -6.32712089705583"
## [1] "2 2 -6.30057977798532"orders <- orders[order(-orders$BIC),]
tail(orders)## m n BIC
## 5 1 0 -6.295251
## 10 2 2 -6.300580
## 8 2 0 -6.308574
## 9 2 1 -6.327121
## 7 1 2 -6.327179
## 6 1 1 -6.353721#1,1
#Final Arima-Garch Model
spec = ugarchspec(variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(3, 0),
include.mean=T), distribution.model="std")
final.model = garchFit(~ arma(3,0)+ garch(1,1), data=gnpgr.train, trace = FALSE)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(final.model)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(3, 0) + garch(1, 1), data = gnpgr.train,
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(3, 0) + garch(1, 1)
## <environment: 0x0000000028d2b128>
## [data = gnpgr.train]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ar2 ar3 omega alpha1
## 9.1582e-03 4.1528e-01 1.8869e-01 -1.5877e-01 1.2275e-05 2.6547e-01
## beta1
## 6.4243e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 9.158e-03 1.527e-03 5.998 2.00e-09 ***
## ar1 4.153e-01 8.377e-02 4.957 7.14e-07 ***
## ar2 1.887e-01 8.075e-02 2.337 0.01946 *
## ar3 -1.588e-01 7.365e-02 -2.156 0.03111 *
## omega 1.228e-05 1.387e-05 0.885 0.37613
## alpha1 2.655e-01 1.077e-01 2.464 0.01373 *
## beta1 6.424e-01 1.973e-01 3.256 0.00113 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 646.3807 normalized: 3.248144
##
## Description:
## Sun Apr 25 15:24:58 2021 by user: perei
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 100.0791 0
## Shapiro-Wilk Test R W 0.9636211 5.207847e-05
## Ljung-Box Test R Q(10) 8.25688 0.603761
## Ljung-Box Test R Q(15) 12.64756 0.6294994
## Ljung-Box Test R Q(20) 25.8959 0.1692776
## Ljung-Box Test R^2 Q(10) 3.452779 0.9686775
## Ljung-Box Test R^2 Q(15) 5.349665 0.988729
## Ljung-Box Test R^2 Q(20) 6.498403 0.9980348
## LM Arch Test R TR^2 4.853983 0.9626789
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.425936 -6.310091 -6.428301 -6.379051#Q3
#Predictions
## Prediction of the return time series
nfore = length(gnpgr.test)
fore.series = NULL
for(f in 1: nfore){
## Fit models
data = gnpgr.train
if(f>=2)
data = c(gnpgr.train,gnpgr.test[1:(f-1)])
final.model = ugarchfit(spec, data, solver = 'hybrid')
## Forecast
fore = ugarchforecast(final.model, n.ahead=1)
fore.series = c(fore.series, fore@forecast$seriesFor)
}
## Compute Accuracy Measures
### Mean Absolute Prediction Error (MAPE)
100*mean(abs(fore.series - gnpgr.test)/abs(gnpgr.test))## [1] 162.224### Precision Measure (PM)
sum((fore.series - gnpgr.test)^2)/sum((gnpgr.test-mean(gnpgr.test))^2)## [1] 0.9352413par(mfrow=c(1,1))
ts.plot(gnpgr.test)
lines(fore.series,col="blue")
Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.