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:

## Libraries
library(TSA)
## Warning: package 'TSA' was built under R version 4.0.4
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(aod)
## Warning: package 'aod' was built under R version 4.0.4
library(urca)
library(tseries)
## Warning: package 'tseries' was built under R version 4.0.4
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(data.table)
## Warning: package 'data.table' was built under R version 4.0.3
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.0.3
library(mgcv)
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
## 
## Attaching package: 'mgcv'
## The following object is masked from 'package:aod':
## 
##     negbin
library(dynlm)
## Warning: package 'dynlm' was built under R version 4.0.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(fpp2)
## Warning: package 'fpp2' was built under R version 4.0.3
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
## -- Attaching packages ------------------------------------------------------------------------- fpp2 2.4 --
## v ggplot2   3.3.0     v fma       2.4  
## v forecast  8.13      v expsmooth 2.3
## Warning: package 'forecast' was built under R version 4.0.3
## Warning: package 'fma' was built under R version 4.0.3
## Warning: package 'expsmooth' was built under R version 4.0.3
## -- Conflicts ---------------------------------------------------------------------------- fpp2_conflicts --
## x forecast::getResponse() masks nlme::getResponse()
library(writexl)
## Warning: package 'writexl' was built under R version 4.0.3
library(ggplot2)
library(feasts)
## Warning: package 'feasts' was built under R version 4.0.3
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.0.3
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
## 
## Attaching package: 'feasts'
## The following object is masked from 'package:nlme':
## 
##     ACF
library(forecast)
library(FitAR)
## Warning: package 'FitAR' was built under R version 4.0.3
## Loading required package: lattice
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 4.0.3
## Loading required package: ltsa
## Warning: package 'ltsa' was built under R version 4.0.3
## Loading required package: bestglm
## Warning: package 'bestglm' was built under R version 4.0.3
## 
## Attaching package: 'FitAR'
## The following object is masked from 'package:forecast':
## 
##     BoxCox
library(lubridate)
## Warning: package 'lubridate' was built under R version 4.0.3
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.0.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
rm(list = ls())
par(mfrow =c(1,1))

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)  

data=data[,4]
data.ts = ts(data,start=c(2010,1,1),freq=365.25)
data.diff = diff(data.ts)
data.diff2 = diff(data.diff)

#Q1 Stationarity Analysis
par(mfrow=c(1,1))
ts.plot(data.ts,col='blue')

acf(data.ts)

ts.plot(data.diff,col='blue')

acf(data.diff)

ts.plot(data.diff2,col='blue')

acf(data.diff2)

#Q2
par(mfrow=c(1,1))
pacf(data.diff)

#ARIMA 1,1,0?

#Q3
#ARIMA order
test_modelA <- function(p,d,q){
  mod = arima(data.ts, 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:3){
  for (d in 0:2){
    for (q in 0:3) {
      possibleError <- tryCatch(
        orders<-rbind(orders,test_modelA(p,d,q)),
        error=function(e) e
      )
      if(inherits(possibleError, "error")) next
      
    }
  }
}
## [1] "0 0 0 19489.8144668527"
## [1] "0 0 1 18527.3691428114"
## [1] "0 0 2 18138.4818734119"
## [1] "0 0 3 17875.7013894288"
## [1] "0 1 0 17836.5479041817"
## [1] "0 1 1 17233.3642636642"
## [1] "0 1 2 17192.3535038164"
## [1] "0 1 3 17190.9309082221"
## [1] "0 2 0 19837.6649459542"
## [1] "0 2 1 17838.1093550864"
## [1] "0 2 2 17237.7931949261"
## [1] "0 2 3 17197.0362338509"
## [1] "1 0 0 17622.0239616309"
## [1] "1 0 1 17229.5820542636"
## [1] "1 0 2 17192.5134600455"
## [1] "1 0 3 17191.7514558669"
## [1] "1 1 0 17524.4697423518"
## [1] "1 1 1 17190.2664706358"
## [1] "1 1 2 17192.246673952"
## [1] "1 1 3 17192.9014158076"
## [1] "1 2 0 18879.5173536383"
## [1] "1 2 1 17526.8384052495"
## [1] "1 2 2 17229.150714451"
## [1] "1 2 3 17197.8289284614"
## [1] "2 0 0 17426.9167624522"
## [1] "2 0 1 17190.9522564778"
## [1] "2 0 2 17192.8912475793"
## [1] "2 0 3 17195.4025477232"
## [1] "2 1 0 17418.998160941"
## [1] "2 1 1 17192.2307793518"
## [1] "2 1 2 17194.2191378882"
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
## [1] "2 1 3 17059.0647928399"
## [1] "2 2 0 18447.3303887929"
## [1] "2 2 1 17421.8334574234"
## [1] "2 2 2 17196.9941616698"
## [1] "2 2 3 17198.9720295063"
## [1] "3 0 0 17359.5293168417"
## [1] "3 0 1 17192.8469038889"
## Warning in log(s2): NaNs produced
## Warning in log(s2): NaNs produced
## [1] "3 0 2 17191.6865417643"
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
## [1] "3 0 3 17182.979276359"
## [1] "3 1 0 17363.4812314533"
## [1] "3 1 1 17186.3821749612"
## [1] "3 1 2 17170.6104601947"
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
## [1] "3 1 3 17053.936918225"
## [1] "3 2 0 18263.7362244792"
## [1] "3 2 1 17366.6567428264"
## [1] "3 2 2 17191.0695029113"
## [1] "3 2 3 17200.9666358973"
orders <- orders[order(-orders$AIC),]
tail(orders)
##    p d q      AIC
## 19 1 1 1 17190.27
## 43 3 1 1 17186.38
## 41 3 0 3 17182.98
## 44 3 1 2 17170.61
## 33 2 1 3 17059.06
## 45 3 1 3 17053.94
arima_final <- arima(data.ts, order=c(3,1,3), method="ML")
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
arima_final
## 
## Call:
## arima(x = data.ts, order = c(3, 1, 3), method = "ML")
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2      ma3
##       1.3448  -1.1203  0.0976  -1.9977  1.9212  -0.7430
## s.e.  0.0344   0.0434  0.0339   0.0255  0.0407   0.0291
## 
## sigma^2 estimated as 325:  log likelihood = -8519.97,  aic = 17051.94
autoplot(arima_final)

#Q4
residuals <- resid(arima_final)
par(mfrow=c(3,1))
ts.plot(residuals)
acf(residuals)
pacf(residuals)

Box.test(resid(arima_final), lag = 7, type = "Box-Pierce", fitdf = 6)
## 
##  Box-Pierce test
## 
## data:  resid(arima_final)
## X-squared = 10.933, df = 1, p-value = 0.0009445
Box.test(resid(arima_final), lag = 7, type = "Ljung-Box", fitdf = 6)
## 
##  Box-Ljung test
## 
## data:  resid(arima_final)
## X-squared = 10.979, df = 1, p-value = 0.0009215
#Q5
vol=time(data.ts)
n = length(data.ts)
nfit = n-8

arima_final = arima(data.ts[0:(nfit)], order=c(3,1,3), method="ML")
## Warning in stats::arima(x = x, order = order, seasonal = seasonal, xreg =
## xreg, : possible convergence problem: optim gave code = 1
arima_pred = as.vector(predict(arima_final,n.ahead=8))
ubound = arima_pred$pred+1.96*arima_pred$se
lbound = arima_pred$pred-1.96*arima_pred$se
ymin = min(lbound)
ymax = max(ubound)

par(mfrow=c(1,1))
plot(vol[nfit:n],data.ts[nfit:n],type="l", ylim=c(ymin,ymax), xlab="Time", main="ARIMA Predictions")
points(vol[(nfit+1):n],arima_pred$pred,col="red")
lines(vol[(nfit+1):n],ubound,lty=3,lwd= 2, col="blue")
lines(vol[(nfit+1):n],lbound,lty=3,lwd= 2, 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.