This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

####Time Series Assignment #####
rm(list = ls())
mydata = read.csv(file.choose(), header = T)
attach(mydata)
library(tseries)
## Warning: package 'tseries' was built under R version 3.4.3
library(TSA)
## Warning: package 'TSA' was built under R version 3.4.3
## Loading required package: leaps
## Loading required package: locfit
## Warning: package 'locfit' was built under R version 3.4.3
## locfit 1.5-9.1    2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 3.4.3
## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(forecast)
## Warning: package 'forecast' was built under R version 3.4.3
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
## 
##     getResponse
library(fpp)                   
## Warning: package 'fpp' was built under R version 3.4.3
## Loading required package: fma
## Warning: package 'fma' was built under R version 3.4.3
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 3.4.3
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.4.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(astsa)
## Warning: package 'astsa' was built under R version 3.4.3
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:fpp':
## 
##     oil
## The following objects are masked from 'package:fma':
## 
##     chicken, sales
## The following object is masked from 'package:forecast':
## 
##     gas
library(DT)
## Warning: package 'DT' was built under R version 3.4.3
names(mydata)
## [1] "OzBeer"
#Plotting a time series for quarterly data. Frequency given in question is give as 4 
beer.ts = ts(OzBeer, frequency = 4)
plot(beer.ts, col= "blue", main= "Quarterly Beer Sales Data", xlab = "Years", ylab= "ML")

#Just looking how the seasons are distributed
seasonplot(beer.ts, year.labels = F, year.labels.left = F, col=1:40, pch=16,
           main= "Quarterly Beer Sales ",
           xlab = "quarters")

#There is some seasonality involved
#Visual representation via box plot
boxplot(beer.ts~cycle(beer.ts), col="blue", main="Quarterly Beer Sales")

p=periodogram(OzBeer)

boxplot(beer.ts~ cycle(beer.ts), main= "Box Plot across time")

#checking if the seasonality is 4
#visual effective using decompose
decompose_beer = decompose(beer.ts, "multiplicative")
plot(as.ts(decompose_beer$seasonal))

plot(as.ts(decompose_beer$trend))

plot(as.ts(decompose_beer$random))

plot(decompose_beer)

#######
#Splitting the screen into 2*2 matrix
par(mfrow = c(2,2))
#Using Moving averages doing EDA
plot(beer.ts, col="gray", main = "1 Year Moving Average Smoothing")
lines(ma(beer.ts, order = 3), col = "red", lwd=3)
plot(beer.ts, col="gray", main = "3 Year Moving Average Smoothing")
lines(ma(beer.ts, order = 12), col = "blue", lwd=3)
plot(beer.ts, col="gray", main = "5 Year Moving Average Smoothing")
lines(ma(beer.ts, order = 36), col = "green", lwd=3)


####Holt Winters Method####
#Doing a as seasonal variation is varying with time
fit = hw(beer.ts, seasonal = "multiplicative")
lines(fitted(fit),col="red")
fitted(fit)
##        Qtr1     Qtr2     Qtr3     Qtr4
## 1  262.8866 226.7882 238.4577 301.2943
## 2  263.5998 225.2712 237.7963 301.5312
## 3  265.1197 227.8819 241.6990 307.3475
## 4  269.9930 231.1087 243.9793 311.0253
## 5  272.8780 235.3873 248.6318 318.2765
## 6  279.0451 241.2209 255.2671 325.9709
## 7  287.6819 247.1007 262.5403 336.2945
## 8  296.9251 255.8731 271.9620 348.5675
## 9  308.5445 267.0340 285.3296 368.1862
## 10 326.1079 283.0524 303.1741 390.5147
## 11 345.6767 298.7547 317.8096 406.3450
## 12 358.0190 308.4920 328.9435 420.7172
## 13 369.3112 320.2329 340.8569 434.7004
## 14 383.3970 330.5917 351.6527 451.1784
## 15 397.3572 341.8388 364.9801 468.5601
## 16 413.3036 356.6229 381.1280 488.4571
## 17 431.0477 371.3347 395.7600 506.5171
## 18 446.5091 386.3690 411.6686 529.5638
#plotting the model
forecast(fit,8)
##       Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 19 Q1       467.4689 449.6426 485.2953 440.2059 494.7320
## 19 Q2       404.1306 388.6368 419.6244 380.4349 427.8264
## 19 Q3       430.8064 414.1248 447.4881 405.2941 456.3188
## 19 Q4       552.4171 530.6915 574.1427 519.1906 585.6435
## 20 Q1       487.6920 468.0895 507.2945 457.7126 517.6715
## 20 Q2       421.4266 404.0029 438.8502 394.7794 448.0738
## 20 Q3       449.0489 429.8322 468.2656 419.6595 478.4384
## 20 Q4       575.5641 549.9233 601.2049 536.3498 614.7783
plot(fit)

fit$model
## Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = beer.ts, seasonal = "multiplicative") 
## 
##   Smoothing parameters:
##     alpha = 0.0643 
##     beta  = 0.0406 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 256.4184 
##     b = 0.022 
##     s=1.1734 0.9247 0.8768 1.0251
## 
##   sigma:  0.0298
## 
##      AIC     AICc      BIC 
## 650.6460 653.5492 671.1360
attributes(fit)
## $names
##  [1] "model"     "mean"      "level"     "x"         "upper"    
##  [6] "lower"     "fitted"    "method"    "series"    "residuals"
## 
## $class
## [1] "forecast"
#AIc - 650.64
#Resetting to a 1*1 matrix
par(mfrow=c(1,2))
#######################

###ARIMA###
View(beer.ts)
kpss.test(beer.ts)
## Warning in kpss.test(beer.ts): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  beer.ts
## KPSS Level = 3.0456, Truncation lag parameter = 1, p-value = 0.01
ndiffs(beer.ts)
## [1] 1
kpss.test(diff(beer.ts, differences = 1))
## Warning in kpss.test(diff(beer.ts, differences = 1)): p-value greater than
## printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff(beer.ts, differences = 1)
## KPSS Level = 0.036853, Truncation lag parameter = 1, p-value = 0.1
#Visual representation of log plot to see how correlation is there in dataset
lag.plot(beer.ts, lags = 16,do.lines = T)

#We can see the correlation is declining from lag 1 to lag 16 
#which makes sense since the correlation normally diminished by time. 
#Seasonal effect is also visible in the plots
#Differencing is done to make the dataset linear
d_beer = diff(log(beer.ts))
#Plotting the differenced dataset
plot(d_beer, main="Differentiated quarterly Sales of Beer")
#Lets do the ACF & PACF plot to identify the "P" & "Q" values
acf2(d_beer)

##         ACF  PACF
##  [1,] -0.09 -0.09
##  [2,] -0.79 -0.81
##  [3,] -0.04 -0.73
##  [4,]  0.88  0.24
##  [5,] -0.06  0.08
##  [6,] -0.76  0.05
##  [7,] -0.03 -0.08
##  [8,]  0.81  0.04
##  [9,] -0.04  0.03
## [10,] -0.73 -0.02
## [11,] -0.01  0.02
## [12,]  0.76  0.10
## [13,] -0.04 -0.01
## [14,] -0.67  0.12
## [15,] -0.03 -0.09
## [16,]  0.73  0.08
## [17,] -0.06 -0.09
## [18,] -0.61  0.09
## [19,] -0.04 -0.04
d.acf = acf(diff(log(beer.ts)))
#Looking at the ACF & PACF Plots, P can be 2 & Q can be 4
#For using in ARIMA, we can use only P or Q
#Lets do some iteration
par(mfrow = c(1,1))

#fitting an arima model
#Iterate arima using q = 0 and p=4
fit1 = arima(log(beer.ts), c(3, 1, 2),seasonal = list(order = c(3, 1, 2), period = 4))
acf(residuals(fit1))

pacf(residuals(fit1))

fit1
## 
## Call:
## arima(x = log(beer.ts), order = c(3, 1, 2), seasonal = list(order = c(3, 1, 
##     2), period = 4))
## 
## Coefficients:
##           ar1      ar2     ar3     ma1      ma2    sar1     sar2     sar3
##       -1.1996  -0.2674  0.0933  0.2438  -0.7561  0.7276  -0.5040  -0.2477
## s.e.   0.1991   0.3207  0.1902  0.1742   0.1632  0.2192   0.1669   0.1840
##          sma1    sma2
##       -1.3408  0.9999
## s.e.   0.2178  0.2116
## 
## sigma^2 estimated as 0.0007511:  log likelihood = 138.42,  aic = -256.84
#AIC = (256.84)
fit2 = arima(log(beer.ts), c(3, 1, 0),seasonal = list(order = c(3, 1, 0), period = 4))
acf(residuals(fit2))

pacf(residuals(fit2))

fit2
## 
## Call:
## arima(x = log(beer.ts), order = c(3, 1, 0), seasonal = list(order = c(3, 1, 
##     0), period = 4))
## 
## Coefficients:
##           ar1      ar2      ar3     sar1     sar2     sar3
##       -1.0607  -0.7264  -0.2884  -0.7469  -0.4792  -0.2845
## s.e.   0.1334   0.1777   0.1748   0.1727   0.1598   0.1394
## 
## sigma^2 estimated as 0.001014:  log likelihood = 133.96,  aic = -255.92
#AIC = (255.92)
#Seeing the AIC chart, the bars are crossing the blue lines
fit3 = arima(log(beer.ts), c(0, 1, 2),seasonal = list(order = c(0, 1, 2), period = 4))
acf(residuals(fit3))

pacf(residuals(fit3))

fit3
## 
## Call:
## arima(x = log(beer.ts), order = c(0, 1, 2), seasonal = list(order = c(0, 1, 
##     2), period = 4))
## 
## Coefficients:
##           ma1     ma2     sma1    sma2
##       -1.0651  0.3171  -0.7356  0.0226
## s.e.   0.1189  0.1220   0.1768  0.1590
## 
## sigma^2 estimated as 0.0009847:  log likelihood = 134.73,  aic = -261.46
#AIC - (261.46)
 
#Conculsion - Fit 1 is the best model
############################
#Some other R functions to forecast Time Series
#using sarima function
par(mfrow=c(1,1))
beer.sarima= sarima(beer.ts,2,1,1,0,1,1,4)
## initial  value 2.856030 
## iter   2 value 2.395107
## iter   3 value 2.336666
## iter   4 value 2.331009
## iter   5 value 2.285166
## iter   6 value 2.277357
## iter   7 value 2.277215
## iter   8 value 2.277023
## iter   9 value 2.276980
## iter  10 value 2.276953
## iter  11 value 2.276952
## iter  12 value 2.276951
## iter  13 value 2.276950
## iter  14 value 2.276950
## iter  15 value 2.276950
## iter  15 value 2.276950
## iter  15 value 2.276950
## final  value 2.276950 
## converged
## initial  value 2.344512 
## iter   2 value 2.343747
## iter   3 value 2.343500
## iter   4 value 2.343452
## iter   5 value 2.343438
## iter   6 value 2.343438
## iter   7 value 2.343437
## iter   7 value 2.343437
## iter   7 value 2.343437
## final  value 2.343437 
## converged

beer.sarima
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1      ar2      ma1     sma1
##       -0.5280  -0.2791  -0.5445  -0.6170
## s.e.   0.1922   0.1700   0.1852   0.1044
## 
## sigma^2 estimated as 103.1:  log likelihood = -252.08,  aic = 514.16
## 
## $degrees_of_freedom
## [1] 63
## 
## $ttable
##      Estimate     SE t.value p.value
## ar1   -0.5280 0.1922 -2.7464  0.0078
## ar2   -0.2791 0.1700 -1.6421  0.1056
## ma1   -0.5445 0.1852 -2.9405  0.0046
## sma1  -0.6170 0.1044 -5.9075  0.0000
## 
## $AIC
## [1] 5.746389
## 
## $AICc
## [1] 5.786793
## 
## $BIC
## [1] 4.87287
beer.sarima$AIC
## [1] 5.746389
#AIC value is 5.74
#Using Auto Arima
arimaauto = auto.arima(beer.ts)
arimaauto
## Series: beer.ts 
## ARIMA(2,1,1)(0,1,1)[4] 
## 
## Coefficients:
##           ar1      ar2      ma1     sma1
##       -0.5280  -0.2791  -0.5445  -0.6170
## s.e.   0.1922   0.1700   0.1852   0.1044
## 
## sigma^2 estimated as 109.6:  log likelihood=-252.08
## AIC=514.16   AICc=515.14   BIC=525.18
#AIC value is 514.16

#Testing Ljung 
Box.test(residuals(fit1), lag = 4, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(fit1)
## X-squared = 0.36498, df = 4, p-value = 0.9852
#Since p value is 0.9188 the data is independantly distributed
#Therefore, we failed to reject the Null hypothesis
#Prediction
arima.predict = predict(fit1, n.ahead = 4*2)
ts.plot(beer.ts,2.718^arima.predict$pred, log="y", lty=c(1,3),
        col="blue",
        main = "Forecast using ARIMA showing Beer Sales",
        xlab="Quarters",
        ylab="ML")

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.