##Basic Set up
library(fpp2)
## Warning: package 'fpp2' was built under R version 4.0.5
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## -- Attaching packages ---------------------------------------------- fpp2 2.4 --
## v ggplot2 3.3.3 v fma 2.4
## v forecast 8.15 v expsmooth 2.3
## Warning: package 'forecast' was built under R version 4.0.5
## Warning: package 'fma' was built under R version 4.0.5
## Warning: package 'expsmooth' was built under R version 4.0.5
##
library(urca)
## Warning: package 'urca' was built under R version 4.0.5
library(fGarch)
## Warning: package 'fGarch' was built under R version 4.0.5
## Loading required package: timeDate
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 4.0.5
## Loading required package: fBasics
## Warning: package 'fBasics' was built under R version 4.0.5
library(vars)
## Warning: package 'vars' was built under R version 4.0.5
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.0.5
##
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
##
## cement, housing, petrol
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.0.5
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.0.5
##
## Attaching package: 'zoo'
## The following object is masked from 'package:timeSeries':
##
## time<-
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 4.0.5
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.0.5
library(hts)
## Warning: package 'hts' was built under R version 4.0.5
setwd("C:\\Users\\samne\\OneDrive\\Desktop\\Master\\School\\Grad\\1st year\\Summer\\Forecasting\\Week 6")
##Download monthly data adjusted close for Microsoft from April 1986 to August 2021
Micro <-read.csv("MSFT-Long.csv", stringsAsFactors= T)
summary(Micro)
## Date Open High Low
## 1986-04-01: 1 Min. : 0.09549 Min. : 0.1068 Min. : 0.09028
## 1986-05-01: 1 1st Qu.: 3.92773 1st Qu.: 4.2168 1st Qu.: 3.82812
## 1986-06-01: 1 Median : 26.57187 Median : 28.0850 Median : 25.12500
## 1986-07-01: 1 Mean : 37.34907 Mean : 39.7656 Mean : 35.55331
## 1986-08-01: 1 3rd Qu.: 37.87500 3rd Qu.: 41.0628 3rd Qu.: 35.49500
## 1986-09-01: 1 Max. :289.48001 Max. :292.9000 Max. :289.31000
## (Other) :420
## Close Adj.Close Volume
## Min. : 0.09809 Min. : 0.06216 Min. :1.825e+07
## 1st Qu.: 4.06445 1st Qu.: 2.57559 1st Qu.:9.194e+08
## Median : 26.61000 Median : 18.54643 Median :1.253e+09
## Mean : 37.98428 Mean : 32.97837 Mean :1.237e+09
## 3rd Qu.: 38.26500 3rd Qu.: 28.42056 3rd Qu.:1.521e+09
## Max. :292.85001 Max. :292.85001 Max. :3.567e+09
##
##Set up a general time series
MSFT = ts(Micro$Adj.Close, start=1986+4/12, frequency = 12)
summary(MSFT)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.06216 2.57559 18.54643 32.97837 28.42056 292.85001
##Plot of Data
autoplot(MSFT) +ggtitle("Microsoft's Montly Adjusted Close from April 1986 to August 2021")
The above plot shows explosive growth in the past decade and small volatility.
##Set up training and validation data
train <- ts(Micro$Adj.Close[1:414], frequency = 12, start = c(1986,4))
test <- ts(Micro$Adj.Close[415:426], frequency = 12, start = c(2020,9))
autoplot(train) + ggtitle("Microsoft's Montly Adjusted Close from April 1986 to August 2020")
autoplot(test) + ggtitle("Microsoft's Montly Adjusted Close from September 2020 to August 2021")
##Seasonality of Data
stl = mstl(train)
autoplot(stl) + ggtitle("STL breakdown of Microsoft Stock")
There is a clear upward trend line and a considerable seasonal pattern that grows in magnitude over time. The remainder also grows considerably, which makes sense as the scale of changes grows greatly as the stock experiences monthly changes that go from less than a tenth of a penny to tens of dollars.
##Auto ETS Model
mod.ets = ets(train)
summary(mod.ets)
## ETS(M,A,N)
##
## Call:
## ets(y = train)
##
## Smoothing parameters:
## alpha = 0.8497
## beta = 0.0216
##
## Initial states:
## l = 0.0593
## b = 0.0102
##
## sigma: 0.0944
##
## AIC AICc BIC
## 2387.579 2387.726 2407.708
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3210141 2.889405 1.476053 0.252459 6.994248 0.2194102
## ACF1
## Training set -0.07971661
##Auto Arima Model
mod.arima = auto.arima(train)
summary(mod.arima)
## Series: train
## ARIMA(1,2,1)
##
## Coefficients:
## ar1 ma1
## -0.3383 -0.9169
## s.e. 0.0526 0.0190
##
## sigma^2 estimated as 7.929: log likelihood=-1011.37
## AIC=2028.75 AICc=2028.81 BIC=2040.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1990218 2.802183 1.47699 0.3450126 7.086425 0.2195495
## ACF1
## Training set 0.002774188
fets = forecast(ets(train), h = 12)
farima = forecast(auto.arima(train), h=12)
autoplot(train) + autolayer(farima, series = "ARIMA Forecast") + autolayer(fets, series = "ETS Forecast", alpha=0.4) + autolayer(test, series = "Test Data", color = "Black", alpha=0.3) +ggtitle("Comparison of Forecasts")
autoplot(test)+ autolayer(fets, series = "ETS Forecast", PI=FALSE) + autolayer(farima, series = "ARIMA Forecast", PI=FALSE)+ggtitle("Comparison of Forecasts to Test Data")
##STL ETS Model
fstle = forecast(stl, method = "ets", h=12)
autoplot(fstle)
##STL Arima
fstla = forecast(stl, method = "arima", h=12)
autoplot(fstla)
autoplot(train) +autolayer(farima, series = "ARIMA Forecast") + autolayer(fets, series = "ETS Forecast", alpha=0.4) + autolayer(fstle, series = "STL ETS Forecast", alpha = 0.3)+ autolayer(fstla, series = "STL ARIMA Forecast", alpha = 0.3)+ autolayer(test, series = "Test Data", color = "Black", alpha=0.3) +ggtitle("Comparison of Forecasts")
autoplot(test) + autolayer(fets, series = "ETS Forecast", PI=FALSE) + autolayer(farima, series = "ARIMA Forecast", PI=FALSE) + autolayer(fstle, series = "STL ETS Forecast", PI=FALSE)+ autolayer(fstla, series = "STL ARIMA Forecast", PI=FALSE) +ggtitle("Comparison of Forecasts to Test Data")
##GARCH models
mod.garch1 <- garchFit(~garch(1,0), data=train)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 0)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 0
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: norm
## h.start: 2
## llh.start: 1
## Length of Series: 414
## Recursion Init: mci
## Series Scale: 36.69453
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -7.28971615 7.289716 0.7289716 TRUE
## omega 0.00000100 100.000000 0.1000000 TRUE
## alpha1 0.00000001 1.000000 0.1000000 TRUE
## gamma1 -0.99999999 1.000000 0.1000000 FALSE
## delta 0.00000000 2.000000 2.0000000 FALSE
## skew 0.10000000 10.000000 1.0000000 FALSE
## shape 1.00000000 10.000000 4.0000000 FALSE
## Index List of Parameters to be Optimized:
## mu omega alpha1
## 1 2 3
## Persistence: 0.1
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 454.96331: 0.728972 0.100000 0.100000
## 1: 392.72699: 0.563942 0.510987 0.983085
## 2: 128.97144: 0.460887 0.0319207 0.960506
## 3: 119.32697: 0.461118 0.0267537 0.960453
## 4: 96.003197: 0.461617 0.0164228 0.960364
## 5: 87.501201: 0.461784 0.0133426 0.960349
## 6: 66.939702: 0.462126 0.00718294 0.960324
## 7: 64.081812: 0.462167 0.00645325 0.960323
## 8: 51.249581: 0.462494 0.000615711 0.960314
## 9: 50.810197: 0.462511 0.00335722 0.960317
## 10: 49.119809: 0.462662 0.00298298 0.961619
## 11: 47.149122: 0.462691 0.00250711 0.961619
## 12: 44.363082: 0.462755 0.00155610 0.961620
## 13: 44.270532: 0.462789 0.00141914 0.961621
## 14: 37.493793: 0.489037 0.00187402 0.962681
## 15: 37.409416: 0.489042 0.00174468 0.962682
## 16: 37.405029: 0.489047 0.00170837 0.962806
## 17: 37.403655: 0.489055 0.00172069 0.962935
## 18: 37.278280: 0.488855 0.00165913 0.998144
## 19: 37.277947: 0.488863 0.00169332 0.998145
## 20: 37.276525: 0.488867 0.00167630 0.998145
## 21: 37.265246: 0.489263 0.00163086 0.998159
## 22: 37.253627: 0.489305 0.00168779 0.997072
## 23: 37.252542: 0.489354 0.00167937 0.995984
## 24: 37.251667: 0.489439 0.00170019 0.993807
## 25: 37.127594: 0.495701 0.00177449 1.00000
## 26: 37.120895: 0.494612 0.00177031 1.00000
## 27: 37.120814: 0.494608 0.00176507 0.999213
## 28: 37.120805: 0.494606 0.00176912 0.998957
## 29: 37.120791: 0.494604 0.00176767 0.999213
## 30: 37.120790: 0.494600 0.00176660 0.999113
## 31: 37.120769: 0.494568 0.00176546 0.999445
## 32: 37.120741: 0.494505 0.00176584 0.999220
## 33: 37.120740: 0.494496 0.00176552 0.999317
## 34: 37.120740: 0.494496 0.00176553 0.999315
##
## Final Estimate of the Negative LLH:
## LLH: 1528.609 norm LLH: 3.692291
## mu omega alpha1
## 18.1453160 2.3772615 0.9993154
##
## R-optimhess Difference Approximated Hessian Matrix:
## mu omega alpha1
## mu -8.010233 3.042810 -1.126269
## omega 3.042810 -36.469018 -6.746849
## alpha1 -1.126269 -6.746849 -158.992549
## attr(,"time")
## Time difference of 0.005985022 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.05086517 secs
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
fgarch1<- predict(mod.garch1, plot=TRUE, n.ahead=12)
fgm1 <- fgarch1$meanForecast
fgm1
## [1] 18.14532 18.14532 18.14532 18.14532 18.14532 18.14532 18.14532 18.14532
## [9] 18.14532 18.14532 18.14532 18.14532
mod.garch2 <- garchFit(~garch(1,1), data=train)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: norm
## h.start: 2
## llh.start: 1
## Length of Series: 414
## Recursion Init: mci
## Series Scale: 36.69453
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -7.28971615 7.289716 0.7289716 TRUE
## omega 0.00000100 100.000000 0.1000000 TRUE
## alpha1 0.00000001 1.000000 0.1000000 TRUE
## gamma1 -0.99999999 1.000000 0.1000000 FALSE
## beta1 0.00000001 1.000000 0.8000000 TRUE
## delta 0.00000000 2.000000 2.0000000 FALSE
## skew 0.10000000 10.000000 1.0000000 FALSE
## shape 1.00000000 10.000000 4.0000000 FALSE
## Index List of Parameters to be Optimized:
## mu omega alpha1 beta1
## 1 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 410.65134: 0.728972 0.100000 0.100000 0.800000
## 1: 342.05890: 0.717856 1.00000e-06 0.107238 0.737826
## 2: 340.31259: 0.717566 0.0625108 0.109015 0.738561
## 3: 225.10021: 0.566462 0.0580810 0.357757 0.467401
## 4: 200.34502: 0.468231 0.0542917 0.284018 0.346154
## 5: 177.55588: 0.418545 1.00000e-06 0.390300 0.270195
## 6: 147.01294: 0.418550 0.0244673 0.390323 0.270217
## 7: 141.53163: 0.491356 0.0245072 0.383656 0.159546
## 8: 113.58320: 0.525883 0.0238252 0.663147 1.00000e-08
## 9: 108.76159: 0.480553 0.0231780 0.751215 1.00000e-08
## 10: 105.33490: 0.494680 0.0224502 0.809893 1.00000e-08
## 11: 100.00870: 0.514895 0.0199635 0.967773 1.00000e-08
## 12: 96.044338: 0.514838 0.0183150 1.00000 1.00000e-08
## 13: 70.700230: 0.512873 0.0104213 1.00000 1.00000e-08
## 14: 66.677878: 0.512868 0.00936015 0.999996 1.00000e-08
## 15: 46.840497: 0.512830 0.000871032 0.999965 1.00000e-08
## 16: 39.695212: 0.512819 0.00280949 0.999966 1.17937e-05
## 17: 39.080014: 0.512531 0.00258523 0.999040 0.00165337
## 18: 38.858021: 0.511406 0.00147775 0.995625 0.00835093
## 19: 37.853579: 0.510269 0.00206464 0.992093 0.0150496
## 20: 33.637334: 0.492442 0.00113033 0.936114 0.123032
## 21: 33.267627: 0.490823 0.000791394 0.919554 0.145351
## 22: 32.381167: 0.492012 0.000763518 0.855414 0.190864
## 23: 31.950898: 0.498671 0.000624356 0.763183 0.250865
## 24: 31.859416: 0.497130 0.000629476 0.775678 0.270150
## 25: 31.823564: 0.495542 0.000663179 0.755455 0.266128
## 26: 31.808206: 0.496261 0.000619023 0.757553 0.270562
## 27: 31.807438: 0.496379 0.000630956 0.757597 0.270343
## 28: 31.807344: 0.496299 0.000627968 0.757039 0.270785
## 29: 31.807343: 0.496308 0.000627744 0.756962 0.270852
## 30: 31.807343: 0.496308 0.000627740 0.756951 0.270866
## 31: 31.807343: 0.496308 0.000627740 0.756949 0.270867
##
## Final Estimate of the Negative LLH:
## LLH: 1523.295 norm LLH: 3.679457
## mu omega alpha1 beta1
## 18.2117775 0.8452442 0.7569489 0.2708674
##
## R-optimhess Difference Approximated Hessian Matrix:
## mu omega alpha1 beta1
## mu -12.186024 -91.31696 1.4582 5.588403
## omega -91.316962 -123.57683 -126.6738 -235.964226
## alpha1 1.458200 -126.67379 -300.5195 -354.924671
## beta1 5.588403 -235.96423 -354.9247 -603.117765
## attr(,"time")
## Time difference of 0.009974003 secs
##
## --- END OF TRACE ---
## Warning in sqrt(diag(fit$cvar)): NaNs produced
##
## Time to Estimate Parameters:
## Time difference of 0.04984212 secs
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
fgarch2<- predict(mod.garch2, plot=TRUE, n.ahead=12)
fgm2 <- fgarch2$meanForecast
fgm2
## [1] 18.21178 18.21178 18.21178 18.21178 18.21178 18.21178 18.21178 18.21178
## [9] 18.21178 18.21178 18.21178 18.21178
These Garch models are largely fruitless as the predictions produced by them are impacted too greatly by the earlier prices of the stock, and thus create poor predictions for the modern upward trend and stock prices. I will not incorporate these two models into the final model evaluation.
##Accuracy Checks and Model Selection
accuracy(fets, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3210141 2.889405 1.476053 0.252459 6.994248 0.2194102
## Test set 23.8063984 30.389392 24.066368 8.747646 8.869968 3.5773835
## ACF1 Theil's U
## Training set -0.07971661 NA
## Test set 0.72686314 2.633552
accuracy(farima, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1990218 2.802183 1.47699 0.3450126 7.086425 0.2195495
## Test set 7.5136609 14.546841 10.74031 2.5029994 3.952733 1.5965110
## ACF1 Theil's U
## Training set 0.002774188 NA
## Test set 0.671489552 1.224348
accuracy(fstle, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3374297 2.590166 1.406848 0.03209148 7.219202 0.2091232
## Test set 21.0252015 26.816578 21.387078 7.72967651 7.899947 3.1791161
## ACF1 Theil's U
## Training set 0.08378838 NA
## Test set 0.69208886 2.340159
accuracy(fstla, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1530267 2.440775 1.376331 0.3946823 7.347004 0.204587
## Test set 1.1737935 7.700812 6.953388 0.2216957 2.711252 1.033597
## ACF1 Theil's U
## Training set 0.004525565 NA
## Test set 0.296937910 0.6773946
The STL ARIMA model had the best performance on the test data across all of the above metrics for error. The ARIMA model also had a strong performance, as it was the second best on all metrics across the four models examined in this section. Ultimately, the STL ARIMA model has the best performance and should be the model selected for use out of the models examined in this assignment to forecast the monthly adjusted close of Microsoft’s stock.