Untitled

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 zoo

library(tseries)

## Warning: package 'tseries' was built under R version 4.0.4

library(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':
## 
##     volatility

library(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':
## 
##     sigma

rm(list = ls())
par(mfrow =c(1,1))


getwd()

## [1] "C:/Users/perei/Dropbox/Worksheets/Lawrence/Python/Time Series Analysis - ISYE-6402-OANO01/Module 3/R Code"

setwd("C:/Users/perei/Dropbox/Worksheets/Lawrence/Python/Time Series Analysis - ISYE-6402-OANO01/Module 3/R Code")
data<-read.csv("TSLA.csv",header=T)


data <- log(data[,2])
data.ts <- ts(data,start=c(2012,1,1),freq=52)
data.growth = diff(data.ts)



#Q1

#Order Selection
test_modelA <- function(p,d,q){
  mod = arima(data.growth, 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:2){
    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 -1001.95153390885"
## [1] "0 0 1 -999.976533799224"
## [1] "0 0 2 -998.781091149895"
## [1] "0 0 3 -997.215363657174"
## [1] "0 0 4 -997.204708031013"
## [1] "0 1 0 -743.007531101551"
## [1] "0 1 1 -993.848790218804"
## [1] "0 1 2 -991.848861232875"
## [1] "0 1 3 -990.468044137523"
## [1] "0 1 4 -988.790283849027"
## [1] "0 2 0 -326.000280490328"
## [1] "0 2 1 -732.121119803844"
## [1] "0 2 2 -974.223227676329"
## [1] "0 2 3 -972.306132578829"
## [1] "0 2 4 -970.722024644068"
## [1] "1 0 0 -999.978980972446"
## [1] "1 0 1 -999.143218168354"
## [1] "1 0 2 -997.816357940132"
## [1] "1 0 3 -995.877880742814"
## [1] "1 0 4 -998.195398558787"
## [1] "1 1 0 -859.69941165643"
## [1] "1 1 1 -991.848866912499"

## Warning in arima(data.growth, order = c(p, d, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## [1] "1 1 2 -990.743965036962"
## [1] "1 1 3 -989.130917058461"
## [1] "1 1 4 -986.482092353511"
## [1] "1 2 0 -556.086756918644"
## [1] "1 2 1 -847.645839092268"
## [1] "1 2 2 -972.253440024993"
## [1] "1 2 3 -971.223059964223"
## [1] "1 2 4 -969.973995873642"
## [1] "2 0 0 -998.913427351276"
## [1] "2 0 1 -997.790707395626"
## [1] "2 0 2 -996.385041232001"
## [1] "2 0 3 -994.40479900847"
## [1] "2 0 4 -1002.0935320065"
## [1] "2 1 0 -901.028953437384"
## [1] "2 1 1 -990.553350539889"
## [1] "2 1 2 -989.130422014283"
## [1] "2 1 3 -987.744041095685"
## [1] "2 1 4 -985.895213595449"
## [1] "2 2 0 -654.99229368305"
## [1] "2 2 1 -888.290921914517"
## [1] "2 2 2 -970.644415138354"

## Warning in log(s2): NaNs produced

## [1] "2 2 3 -969.330067828764"
## [1] "2 2 4 -967.134770596984"
## [1] "3 0 0 -997.256885851907"
## [1] "3 0 1 -995.841274461945"
## [1] "3 0 2 -994.392863657056"

## Warning in log(s2): NaNs produced

## [1] "3 0 3 -992.855109486565"
## [1] "3 0 4 -1000.31583671188"
## [1] "3 1 0 -930.317317023143"
## [1] "3 1 1 -988.795182282954"
## [1] "3 1 2 -987.135492150675"
## [1] "3 1 3 -985.869986436799"

## Warning in arima(data.growth, order = c(p, d, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## [1] "3 1 4 -984.312930853974"
## [1] "3 2 0 -743.08402119346"
## [1] "3 2 1 -917.006672695157"
## [1] "3 2 2 -968.771706666303"
## [1] "3 2 3 -967.70219128848"

## Warning in log(s2): NaNs produced

## [1] "3 2 4 -965.256954878143"
## [1] "4 0 0 -997.275179594166"
## [1] "4 0 1 -996.907630304939"
## [1] "4 0 2 -1003.80278179236"

## Warning in arima(data.growth, order = c(p, d, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## [1] "4 0 3 -1001.77813944674"
## [1] "4 0 4 -998.32022338631"
## [1] "4 1 0 -932.817470576436"
## [1] "4 1 1 -988.643859271108"
## [1] "4 1 2 -988.300379812244"
## [1] "4 1 3 -995.302905244194"

## Warning in arima(data.growth, order = c(p, d, q), method = "ML"): possible
## convergence problem: optim gave code = 1

## [1] "4 1 4 -993.345114751789"
## [1] "4 2 0 -768.824540958537"
## [1] "4 2 1 -919.294422710676"
## [1] "4 2 2 -968.663006972023"

## Warning in log(s2): NaNs produced

## [1] "4 2 3 -968.203943910775"
## [1] "4 2 4 -967.606352267544"

orders <- orders[order(-orders$AIC),]
tail(orders)

##    p d q       AIC
## 17 1 0 0  -999.979
## 51 3 0 4 -1000.316
## 65 4 0 3 -1001.778
## 2  0 0 0 -1001.952
## 36 2 0 4 -1002.094
## 64 4 0 2 -1003.803

#4,0,2

model.arima = arima(data.growth, order=c(4,0,2), method="ML")
model.arima

## 
## Call:
## arima(x = data.growth, order = c(4, 0, 2), method = "ML")
## 
## Coefficients:
##           ar1      ar2     ar3     ar4     ma1     ma2  intercept
##       -1.2930  -0.8478  0.1053  0.1200  1.3200  0.9343     0.0058
## s.e.   0.0662   0.1212  0.0848  0.0573  0.0446  0.0928     0.0036
## 
## sigma^2 estimated as 0.00405:  log likelihood = 509.9,  aic = -1003.8

#Residuals
resids = resid(model.arima)

#Residual ACF PLots
par(mfrow=c(1,2))
acf(resids,main='Residual ACF')
acf(resids^2,main='Squared Residual ACF')

# for serial correlation
Box.test(resids,lag=7,type='Ljung',fitdf=6)

## 
##  Box-Ljung test
## 
## data:  resids
## X-squared = 0.883, df = 1, p-value = 0.3474

Box.test((resids)^2,lag=7,type='Ljung',fitdf=6)

## 
##  Box-Ljung test
## 
## data:  (resids)^2
## X-squared = 6.9138, df = 1, p-value = 0.008553

#Q2

#Initial GARCH Order
#ARIMA-GARCH GARCH order
test_modelAGG <- function(m,n){
  spec = ugarchspec(variance.model=list(garchOrder=c(m,n)),
                    mean.model=list(armaOrder=c(4,2), 
                                    include.mean=T), distribution.model="std")    
  fit = ugarchfit(spec, data.growth, 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:3){
  for (n in 0:3){
    possibleError <- tryCatch(
      orders<-rbind(orders,test_modelAGG(m,n)),
      error=function(e) e
    )
    if(inherits(possibleError, "error")) next
  }
}

## [1] "0 0 -2.56390336438806"
## [1] "0 1 -2.54828047735607"
## [1] "0 2 -2.52640253908174"
## [1] "0 3 -2.51717388475447"
## [1] "1 0 -2.44543643661524"
## [1] "1 1 -2.54946776356549"
## [1] "1 2 -2.58530066319158"
## [1] "1 3 -2.58614144586915"
## [1] "2 0 -2.53658372933802"
## [1] "2 1 -2.53390384643259"
## [1] "2 2 -2.51887509720818"
## [1] "2 3 -2.56877973208704"
## [1] "3 0 -2.50784345663187"
## [1] "3 1 -2.57499071854611"
## [1] "3 2 -2.55572211886978"
## [1] "3 3 -2.55683765724558"

orders <- orders[order(-orders$BIC),]
tail(orders)

##    m n       BIC
## 17 3 3 -2.556838
## 2  0 0 -2.563903
## 13 2 3 -2.568780
## 15 3 1 -2.574991
## 8  1 2 -2.585301
## 9  1 3 -2.586141

#1,2

#ARMA update
#ARIMA-GARCH ARIMA order
test_modelAGA <- function(p,q){
  spec = ugarchspec(variance.model=list(garchOrder=c(1,3)),
                    mean.model=list(armaOrder=c(p,q), 
                                    include.mean=T), distribution.model="std")    
  fit = ugarchfit(spec, data.growth, 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 -2.58397060138191"
## [1] "0 1 -2.5690970307804"
## [1] "0 2 -2.55453618381674"
## [1] "0 3 -2.53914769890282"
## [1] "0 4 -2.53058023137355"
## [1] "1 0 -2.56915169555075"
## [1] "1 1 -2.55726242932836"
## [1] "1 2 -2.54198813263015"
## [1] "1 3 -2.5268975978481"
## [1] "1 4 -2.52426342653738"
## [1] "2 0 -2.55410280972165"
## [1] "2 1 -2.53887935632634"
## [1] "2 2 -2.54925085559696"
## [1] "2 3 -2.5687482770954"
## [1] "2 4 -2.51839919123002"
## [1] "3 0 -2.53975699297666"
## [1] "3 1 -2.52534822284405"
## [1] "3 2 -2.5116127099138"
## [1] "3 3 -2.50947148147104"
## [1] "3 4 -2.50398310989308"
## [1] "4 0 -2.53118923211884"
## [1] "4 1 -2.52039448099236"
## [1] "4 2 -2.58614144586915"
## [1] "4 3 -2.57789023732717"
## [1] "4 4 -2.54510585514929"

orders <- orders[order(-orders$BIC),]
tail(orders)

##    p q       BIC
## 15 2 3 -2.568748
## 3  0 1 -2.569097
## 7  1 0 -2.569152
## 25 4 3 -2.577890
## 2  0 0 -2.583971
## 24 4 2 -2.586141

#4,4

#GARCH update
test_modelAGG <- function(m,n){
  spec = ugarchspec(variance.model=list(garchOrder=c(m,n)),
                    mean.model=list(armaOrder=c(4,2), 
                                    include.mean=T), distribution.model="std")    
  fit = ugarchfit(spec, data.growth, 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:3){
  for (n in 0:3){
    possibleError <- tryCatch(
      orders<-rbind(orders,test_modelAGG(m,n)),
      error=function(e) e
    )
    if(inherits(possibleError, "error")) next
  }
}

## [1] "0 0 -2.56390336438806"
## [1] "0 1 -2.54828047735607"
## [1] "0 2 -2.52640253908174"
## [1] "0 3 -2.51717388475447"
## [1] "1 0 -2.44543643661524"
## [1] "1 1 -2.54946776356549"
## [1] "1 2 -2.58530066319158"
## [1] "1 3 -2.58614144586915"
## [1] "2 0 -2.53658372933802"
## [1] "2 1 -2.53390384643259"
## [1] "2 2 -2.51887509720818"
## [1] "2 3 -2.56877973208704"
## [1] "3 0 -2.50784345663187"
## [1] "3 1 -2.57499071854611"
## [1] "3 2 -2.55572211886978"
## [1] "3 3 -2.55683765724558"

orders <- orders[order(-orders$BIC),]
tail(orders)

##    m n       BIC
## 17 3 3 -2.556838
## 2  0 0 -2.563903
## 13 2 3 -2.568780
## 15 3 1 -2.574991
## 8  1 2 -2.585301
## 9  1 3 -2.586141

#1,1

#Final order (4,4)x(1,1)
final.model = garchFit(~ arma(4,4)+ garch(1,1), data=data.growth, trace = FALSE)

## Warning in arima(.series$x, order = c(u, 0, v), include.mean = include.mean):
## possible convergence problem: optim gave code = 1

## 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(4, 4) + garch(1, 1), data = data.growth, 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(4, 4) + garch(1, 1)
## <environment: 0x0000000027d69d60>
##  [data = data.growth]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu          ar1          ar2          ar3          ar4          ma1  
##  0.00426134  -0.75959668  -0.19646790   0.57980707   0.09368805   0.76189873  
##         ma2          ma3          ma4        omega       alpha1        beta1  
##  0.20718093  -0.54108488   0.01292577   0.00026663   0.06429873   0.87028504  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      0.0042613   0.0049026    0.869  0.38474    
## ar1    -0.7595967   0.0848621   -8.951  < 2e-16 ***
## ar2    -0.1964679   0.1382944   -1.421  0.15542    
## ar3     0.5798071   0.1484487    3.906 9.39e-05 ***
## ar4     0.0936881   0.1006034    0.931  0.35172    
## ma1     0.7618987   0.1004590    7.584 3.35e-14 ***
## ma2     0.2071809   0.1584895    1.307  0.19114    
## ma3    -0.5410849   0.1672439   -3.235  0.00122 ** 
## ma4     0.0129258   0.1073348    0.120  0.90415    
## omega   0.0002666   0.0001437    1.856  0.06352 .  
## alpha1  0.0642987   0.0311171    2.066  0.03880 *  
## beta1   0.8702850   0.0469443   18.539  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  522.7252    normalized:  1.368391 
## 
## Description:
##  Sun Apr 25 13:17:59 2021 by user: perei 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  55.67742  8.124612e-13
##  Shapiro-Wilk Test  R    W      0.9822455 0.0001214266
##  Ljung-Box Test     R    Q(10)  4.683363  0.9113019   
##  Ljung-Box Test     R    Q(15)  5.972257  0.98021     
##  Ljung-Box Test     R    Q(20)  8.220803  0.9903032   
##  Ljung-Box Test     R^2  Q(10)  2.410858  0.992112    
##  Ljung-Box Test     R^2  Q(15)  5.080218  0.9914252   
##  Ljung-Box Test     R^2  Q(20)  6.10085   0.9987543   
##  LM Arch Test       R    TR^2   4.568772  0.9708563   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -2.673954 -2.550014 -2.675849 -2.624784

#Q3
train = data.growth[1:(length(data.growth)-40)]
test = data.growth[(length(data.growth)-39):length(data.growth)]

nfore = length(test)
fore.series = NULL

spec = ugarchspec(variance.model=list(garchOrder=c(1,1)),
                    mean.model=list(armaOrder=c(4, 4), 
                                    include.mean=T), distribution.model="std")    
final.model = garchFit(~ arma(4,4)+ garch(1,1), data=train, trace = FALSE)

## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.

for(f in 1: nfore){
  ## Fit models
  data = train
  if(f>=2)
    data = c(train,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)
  
}

#Plots
par(mfrow=c(1,1))
ts.plot(test,main='TSLA Growth Forecasts',ylab='Growth')
lines(fore.series,col="blue")

#Accuracy Measures
preds = as.vector(fore.series)
obs = as.vector(test)


#Mean Absolute Percentage Measure (MAPE)
100*mean(abs(preds-obs)/abs(obs))

## [1] 311.6466

#311.66

#Precision
sum((preds-obs)^2)/sum((obs-mean(obs))^2)

## [1] 1.136645

#1.14

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.