## Libraries
library(fastDummies)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(generics)
##
## Attaching package: 'generics'
## The following object is masked from 'package:caret':
##
## train
## The following object is masked from 'package:lubridate':
##
## as.difftime
## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, union
library(tsibble)
##
## Attaching package: 'tsibble'
## The following object is masked from 'package:lubridate':
##
## interval
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:generics':
##
## explain
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(fpp3)
## ── Attaching packages ──────────────────────────────────────────── fpp3 0.4.0 ──
## ✓ tibble 3.1.6 ✓ feasts 0.2.2
## ✓ tidyr 1.2.0 ✓ fable 0.3.1
## ✓ tsibbledata 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks generics::intersect(), base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x fabletools::MAE() masks caret::MAE()
## x fabletools::RMSE() masks caret::RMSE()
## x tsibble::setdiff() masks generics::setdiff(), base::setdiff()
## x tsibble::union() masks generics::union(), base::union()
library(modeest)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Registered S3 method overwritten by 'forecast':
## method from
## predict.default statip
##
## Attaching package: 'forecast'
## The following object is masked from 'package:modeest':
##
## naive
## The following objects are masked from 'package:generics':
##
## accuracy, forecast
library(latex2exp)
library(seasonal)
##
## Attaching package: 'seasonal'
## The following object is masked from 'package:tibble':
##
## view
library(rugarch)
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:fabletools':
##
## report
## The following object is masked from 'package:stats':
##
## sigma
## Data Set
df <- read.csv("//Users//kevinclifford//Downloads//Alcohol_Sales.csv", header=TRUE)
df$Sales <- df$S4248SM144NCEN
df$S4248SM144NCEN <- NULL
ts <- ts(df$Sales, frequency = 12, start=c(1992))
plot(ts)
## ETS Models
fit1 <- ets(ts)
fit1
## ETS(M,Ad,M)
##
## Call:
## ets(y = ts)
##
## Smoothing parameters:
## alpha = 0.0805
## beta = 0.0232
## gamma = 1e-04
## phi = 0.9592
##
## Initial states:
## l = 4199.083
## b = 3.4466
## s = 1.1642 1.0362 1.0338 0.9829 1.0534 1.0081
## 1.106 1.0665 0.9754 0.9758 0.8275 0.7702
##
## sigma: 0.0455
##
## AIC AICc BIC
## 5672.314 5674.549 5740.422
plot(fit1)
accuracy(fit1)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 53.45892 364.3906 285.5164 0.5004852 3.657119 0.6753521 -0.2918586
fe <- forecast(fit1, 12)
acc <- accuracy(fe, df$Sales[1:12])
acc
## ME RMSE MAE MPE MAPE MASE
## Training set 53.45892 364.3906 285.5164 0.5004852 3.657119 0.312568
## Test set -9484.70207 9600.8980 9484.7021 -228.2510458 228.251046 10.383342
## ACF1
## Training set -0.2918586
## Test set NA
plot(fe, main="MMN")
arima1 <- auto.arima(ts)
arima1
## Series: ts
## ARIMA(3,1,1)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 ma1 sma1 sma2
## -0.1428 0.1580 0.5125 -0.9483 -0.2601 -0.2642
## s.e. 0.0637 0.0651 0.0609 0.0328 0.0581 0.0543
##
## sigma^2 = 102379: log likelihood = -2242.28
## AIC=4498.56 AICc=4498.93 BIC=4524.77
plot(arima1)
accuracy(arima1)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 34.22564 310.4741 232.8522 0.3437269 2.939946 0.5507817 0.02751723
fe2 <- forecast(arima1, 12)
acc2 <- accuracy(fe2, df$Sales[1:12])
acc2
## ME RMSE MAE MPE MAPE MASE
## Training set 34.22564 310.4741 232.8522 0.3437269 2.939946 0.2549141
## Test set -9644.03404 9745.1509 9644.0340 -232.2911945 232.291195 10.5577699
## ACF1
## Training set 0.02751723
## Test set NA
plot(fe2, main="Auto-ARIMA")
garch_spec <- ugarchspec(mean.model=list(armaOrder = c(1,0)))
garch_train<- ugarchfit(garch_spec, data =ts)
plot(garch_train, which = "all")
##
## please wait...calculating quantiles...
garch_train
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 5.3102e+03 6.2360e+02 8.51544 0.000000
## ar1 9.0107e-01 2.7546e-02 32.71121 0.000000
## omega 8.5115e+03 9.6655e+03 0.88061 0.378530
## alpha1 5.1884e-02 1.4424e-02 3.59696 0.000322
## beta1 9.4712e-01 1.9085e-02 49.62720 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 5.3102e+03 6.4039e+02 8.2921 0.000000
## ar1 9.0107e-01 2.3223e-02 38.7998 0.000000
## omega 8.5115e+03 8.4620e+03 1.0059 0.314486
## alpha1 5.1884e-02 1.2988e-02 3.9947 0.000065
## beta1 9.4712e-01 1.7599e-02 53.8159 0.000000
##
## LogLikelihood : -2759.017
##
## Information Criteria
## ------------------------------------
##
## Akaike 17.009
## Bayes 17.068
## Shibata 17.009
## Hannan-Quinn 17.033
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 43.50 4.24e-11
## Lag[2*(p+q)+(p+q)-1][2] 44.06 0.00e+00
## Lag[4*(p+q)+(p+q)-1][5] 47.59 0.00e+00
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 3.432 0.063942
## Lag[2*(p+q)+(p+q)-1][5] 7.996 0.029537
## Lag[4*(p+q)+(p+q)-1][9] 14.911 0.003759
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 4.155 0.500 2.000 0.041514
## ARCH Lag[5] 8.674 1.440 1.667 0.013805
## ARCH Lag[7] 12.526 2.315 1.543 0.004638
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.2275
## Individual Statistics:
## mu 0.8109
## ar1 0.1799
## omega 0.1555
## alpha1 0.1453
## beta1 0.2785
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.9564 3.396e-01
## Negative Sign Bias 0.6179 5.371e-01
## Positive Sign Bias 4.5384 8.032e-06 ***
## Joint Effect 25.9659 9.696e-06 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 104.2 9.113e-14
## 2 30 126.3 4.163e-14
## 3 40 137.4 6.887e-13
## 4 50 158.2 1.992e-13
##
##
## Elapsed time : 0.1139069
forecast_garch <- ugarchforecast(garch_train, n.ahead = 4)
plot(forecast_garch, which = 3)
