##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.