Project #2
Charles Bickel (305012336)
6/3/2020
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls(all=TRUE))
# Load Libraries
library(lattice)
library(foreign)
library(MASS)
library(car)## Loading required package: carDatarequire(stats)
require(stats4)## Loading required package: stats4library(KernSmooth)## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009library(fastICA)
library(cluster)
library(leaps)
library(mgcv)## Loading required package: nlme## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.library(rpart)
library(pan)
library(mgcv)
library(DAAG)##
## Attaching package: 'DAAG'## The following object is masked from 'package:car':
##
## vif## The following object is masked from 'package:MASS':
##
## hillslibrary("TTR")
library(tis)##
## Attaching package: 'tis'## The following object is masked from 'package:TTR':
##
## lags## The following object is masked from 'package:mgcv':
##
## tirequire("datasets")
require(graphics)
library("forecast")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'forecast'## The following object is masked from 'package:tis':
##
## easter## The following object is masked from 'package:nlme':
##
## getResponse#install.packages("astsa")
#require(astsa)
library(xtable)
# New libraries added:
library(TSA)## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(timeSeries)## Loading required package: timeDate##
## Attaching package: 'timeDate'## The following objects are masked from 'package:TSA':
##
## kurtosis, skewness## The following object is masked from 'package:xtable':
##
## align## The following objects are masked from 'package:tis':
##
## dayOfWeek, dayOfYear, isHoliday##
## Attaching package: 'timeSeries'## The following objects are masked from 'package:tis':
##
## description, interpNAlibrary(fUnitRoots)## Loading required package: fBasics##
## Attaching package: 'fBasics'## The following object is masked from 'package:TTR':
##
## volatility## The following object is masked from 'package:car':
##
## densityPlotlibrary(fBasics)
library(tseries)
library(timsac)
library(TTR)
library(fpp)## Loading required package: fma##
## Attaching package: 'fma'## The following objects are masked from 'package:DAAG':
##
## milk, ozone## The following objects are masked from 'package:MASS':
##
## cement, housing, petrol## Loading required package: expsmooth## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(vars)## Loading required package: strucchange## Loading required package: sandwich## Loading required package: urca##
## Attaching package: 'urca'## The following objects are masked from 'package:fUnitRoots':
##
## punitroot, qunitroot, unitrootTablelibrary(lmtest)
library(dlnm)## This is dlnm 2.3.9. For details: help(dlnm) and vignette('dlnmOverview').library(strucchange)
library(ggplot2)
library(seasonal)##
## Attaching package: 'seasonal'## The following objects are masked from 'package:timeSeries':
##
## outlier, seriesPart (a)
The data I will be looking at and modeling are Apple and Microsoft daily stock prices from March of 1986. I expect the time series to share similar dynamics because they are both tech giants in the United States market and compete in the same market. Both time series grow with an exponential or quadratic trend and begin to grow very steeply around 2016/17. We can see that both companies stock prices take a drop around the 2009 financial crisis. Both series show strong signals of seasonality and cycles. Both time series do not look to be stationary as there is very little mean reversion. The ACF and PACF plots are messy to read because of our high order of frequency, but we can tell that the data does not look covariance stationary as the ACF decays very slowly from around 1 in each time series. The PACF plots look like they may resembe an AR or S-AR process so an ARIMA model may be good fit for the data. Before beginning modeling, we take the first difference of both time series to stabalize the volatilities and make the data covariance stationary. Before taking the log, I squared the data so that the values are not so small that we are looking at Looking at the time series plots of the transformed data, there is strong mean reversion which suggests our data is now covariace stationary. The volatility has also been tamed quiet a bit as the range of values has dropped by a factor of 100. Now looking at the ACF and PACF plots, it is hard to tell what process we are looking it but it looks similiar to an S-ARIMA model with clear seasonal patterns in the lag spikes. The ACF decays much faster to 0 now, and on the contrary, the PACFs decay much slower. Interesting note, we can see that Microsoft had a higher stock price than Apple until 2009, when the financial crisis hit. Our first differenced data appears to have multiplicative seasonality.
#load in data
msft <- read.csv("Data/MSFT_Stock.csv")
appl <- read.csv("Data/Apple_stock.csv")
head(msft)## MacroTrends.Data.Download
## 1 MSFT - Historical Price and Volume Data
## 2 Note: Historical prices are adjusted for stock splits.
## 3
## 4 Disclaimer and Terms of Use: Historical stock data is provided 'as is' and solely for informational purposes, not for trading purposes or advice.
## 5 MacroTrends LLC expressly disclaims the accuracy, adequacy, or completeness of any data and shall not be liable for any errors, omissions or other defects in,
## 6 delays or interruptions in such data, or for any actions taken in reliance thereon. Neither MacroTrends LLC nor any of our information providers will be liable
## X X.1 X.2 X.3 X.4
## 1
## 2
## 3
## 4
## 5
## 6head(appl)## MacroTrends.Data.Download
## 1 AAPL - Historical Price and Volume Data
## 2 Note: Historical prices are adjusted for stock splits.
## 3
## 4 Disclaimer and Terms of Use: Historical stock data is provided 'as is' and solely for informational purposes, not for trading purposes or advice.
## 5 MacroTrends LLC expressly disclaims the accuracy, adequacy, or completeness of any data and shall not be liable for any errors, omissions or other defects in,
## 6 delays or interruptions in such data, or for any actions taken in reliance thereon. Neither MacroTrends LLC nor any of our information providers will be liable
## X X.1 X.2 X.3 X.4
## 1
## 2
## 3
## 4
## 5
## 6tail(msft)## MacroTrends.Data.Download X X.1 X.2 X.3 X.4
## 8636 5/27/2020 180.2 181.9877 176.6 181.81 39517146
## 8637 5/28/2020 180.74 184.1474 180.38 181.4 33526412
## 8638 5/29/2020 182.73 184.27 180.41 183.25 42146720
## 8639 6/1/2020 182.54 183 181.46 182.83 22668821
## 8640 6/2/2020 184.25 185 181.35 184.91 30674990
## 8641 6/3/2020 184.815 185.94 183.58 185.36 26441229tail(appl)## MacroTrends.Data.Download X X.1 X.2 X.3 X.4
## 9961 5/27/2020 316.14 318.71 313.09 318.11 28236274
## 9962 5/28/2020 316.77 323.44 315.63 318.25 33281977
## 9963 5/29/2020 319.25 321.15 316.47 317.94 38399532
## 9964 6/1/2020 317.75 322.35 317.21 321.85 20254653
## 9965 6/2/2020 320.745 323.44 318.93 323.34 21754815
## 9966 6/3/2020 324.66 326.2 322.3 325.12 25574067appl_ts <- ts(as.numeric(appl$X.3[seq(1340, 9966)]), start = c(1986, 49), frequency = 252)
head(appl_ts)## Time Series:
## Start = c(1986, 49)
## End = c(1986, 54)
## Frequency = 252
## [1] 0.4420 0.4664 0.4643 0.4798 0.4732 0.5045msft_ts <- ts(as.numeric(msft$X.3[seq(15, 8641)]), start = c(1986, 49), frequency = 252)
head(msft_ts)## Time Series:
## Start = c(1986, 49)
## End = c(1986, 54)
## Frequency = 252
## [1] 0.0972 0.1007 0.1024 0.0998 0.0981 0.0955length(msft_ts)## [1] 8627length(appl_ts)## [1] 8627#Show plot of time series
plot(appl_ts, main = "Apple Stock Closing Prices(Daily)")
plot(msft_ts, main = "Microsoft Stock Closing Prices(Daily)")
plot(appl_ts, main = "Apple Vs. Microsoft Historical Stock Prices")
nberShade()
grid()
lines(msft_ts, col = "Blue")
lines(appl_ts)
#show acf & pacf plots
par(mfrow = c(2,2))
plot(acf(appl_ts, lag.max = 50))
plot(pacf(appl_ts, lag.max = 50))
plot(acf(msft_ts, lag.max = 50))
plot(pacf(msft_ts, lag.max = 50))
#take first logged difference of our squared data sets so we can get data stationary
APPL <- diff(appl_ts)
MSFT <- diff(msft_ts)
par(mfrow = c(1,1))
autoplot(APPL) + ggtitle("Apple Stock Closing Prices(d=1)")
autoplot(MSFT) + ggtitle("Microsoft Stock Closing Prices(d=1)")
ggtsdisplay(APPL, lag.max = 50) + ggtitle("Apple Stock Closing Prices(d=1)")
ggtsdisplay(MSFT, lag.max = 50) + ggtitle("Microsoft Stock Closing Prices(d=1)")
plot(APPL, main = "Apple Vs. Microsoft Historical Stock Prices")
nberShade()
grid()
lines(MSFT, col = "Blue")
lines(APPL)
plot(APPL, main = "Apple Vs. Microsoft Historical Stock Prices", xlim = c(2009.4,2009.75), ylim = c(-2, 2))
nberShade()
grid()
lines(MSFT, col = "Blue")
lines(APPL)
Part (b)
We are going to create a few different forecasting models for our time series. The first model we are going to create is a model with trend and seasonal components. We start by splitting up our data into a rough 2/3 split into training and testing. We went a little over 2/3 for our training so that way we are able to capture a little more of the behavior once the stock prices began to exponentially rise. Once we create our seasonally adjusted model and forecast, we will use the seasonally adjusted data to then account for cycles using and ARIMA model, ets model, and a Holts-Winters model.
#split data up into training and testing
APPL_training <- window(APPL, start = c(1986, 50), end = c(2014, 168))
APPL_test <- window(APPL, start = c(2014,169))
MSFT_training <- window(MSFT, start = c(1986, 50), end = c(2014, 168))
MSFT_test <- window(MSFT, start = c(2014,169))
#create time series dumby variable
t <- seq(1, length(APPL_training))
#seasonal decompositions
Afit1 <- stl(APPL_training, t.window = 13, s.window = "periodic", robust = TRUE)
autoplot(Afit1) + ggtitle("Stl Decomposition of Apple Stock Prices")
Mfit1 <- stl(MSFT_training, t.window = 13, s.window = "periodic", robust = TRUE)
autoplot(Mfit1) + ggtitle("Stl Decomposition of Microsoft Stock Prices")
summary(Afit1)## Call:
## stl(x = APPL_training, s.window = "periodic", t.window = 13,
## robust = TRUE)
##
## Time.series components:
## seasonal trend remainder
## Min. :-0.05144448 Min. :-1.112892 Min. :-8.725975
## 1st Qu.:-0.01187283 1st Qu.:-0.018840 1st Qu.:-0.030584
## Median : 0.00172237 Median :-0.000064 Median : 0.000000
## Mean : 0.00013557 Mean : 0.005810 Mean : 0.008227
## 3rd Qu.: 0.01278741 3rd Qu.: 0.024670 3rd Qu.: 0.031358
## Max. : 0.04280597 Max. : 3.316474 Max. : 7.354725
## IQR:
## STL.seasonal STL.trend STL.remainder data
## 0.02466 0.04351 0.06194 0.08030
## % 30.7 54.2 77.1 100.0
##
## Weights:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.4397 0.9452 0.7061 0.9910 1.0000
##
## Other components: List of 5
## $ win : Named num [1:3] 71751 13 253
## $ deg : Named int [1:3] 0 1 1
## $ jump : Named num [1:3] 7176 2 26
## $ inner: int 1
## $ outer: int 15summary(Mfit1)## Call:
## stl(x = MSFT_training, s.window = "periodic", t.window = 13,
## robust = TRUE)
##
## Time.series components:
## seasonal trend remainder
## Min. :-0.15164559 Min. :-1.1606713 Min. :-6.771198
## 1st Qu.:-0.02675918 1st Qu.:-0.0522197 1st Qu.:-0.076309
## Median : 0.00157607 Median : 0.0010881 Median : 0.002083
## Mean :-0.00000140 Mean :-0.0056318 Mean : 0.011873
## 3rd Qu.: 0.02715653 3rd Qu.: 0.0474844 3rd Qu.: 0.084560
## Max. : 0.11377104 Max. : 0.7503258 Max. : 4.945617
## IQR:
## STL.seasonal STL.trend STL.remainder data
## 0.05392 0.09970 0.16087 0.23000
## % 23.4 43.3 69.9 100.0
##
## Weights:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.4900 0.9452 0.7116 0.9926 1.0000
##
## Other components: List of 5
## $ win : Named num [1:3] 71751 13 253
## $ deg : Named int [1:3] 0 1 1
## $ jump : Named num [1:3] 7176 2 26
## $ inner: int 1
## $ outer: int 15#forecast seasonally adjusted models
Afit1 %>% forecast(method = "naive") %>%
autoplot() + ggtitle("APPL Stock Forecasts from STL + Random Walk")
Mfit1 %>% forecast(method = "naive") %>%
autoplot() + ggtitle("MSFT Stock Forecasts from STL + Random Walk")
# Afit1 <- tslm(APPL_training ~ trend + season)
# Mfit1 <- tslm(MSFT_training ~ trend + season)I attempted to use auto.arima to calculate the best ARIMA models to use for modeling our seasonally adjusted series. However, when attempting to run auto.arima() on the time series, every time it crashed RStudio. You can see the code that I attempted commented out. Therefore, I had to manually build the ARMA models using arma(). I used autoarmafit() to find the R’s best predicted fit and for APPL Stocks I used an ARMA(5,3) and for MSFT Stocks I used an ARMA(3,2). For APPL, the AR(2), AR(3), AR(5), and MA(2) are significant at a 95% level and should be used in the ARMA(5,3) model, while for MSFT, the AR(3) & MA(1) were significant and should be used in the ARMA(3,2) model. The summary statistics show us getting extremely low AIC and BIC values of around -20,000 which are very good and suggest that these are good forecasting models. Note, even though the MA(2) is not statistically significant in the ARMA(3,2) for MSFT Stocks, the AIC and BIC values were ever so slighlty lower for the ARMA(3,2) than and ARMA(3,1) model.
#build S-ARIMA model
APPL_training %>% stl(s.window = "periodic") %>% seasadj() -> APPL_training_SA
APPL_test %>% stl(s.window = "periodic") %>% seasadj() -> APPL_test_SA
APPL %>% stl(s.window = "periodic") %>% seasadj() -> APPL_SA
autoplot(APPL_training_SA) + ggtitle("Seasonally Adjusted Apple Prices")
MSFT_training %>% stl(s.window = "periodic") %>% seasadj() -> MSFT_training_SA
MSFT_test %>% stl(s.window = "periodic") %>% seasadj() -> MSFT_test_SA
MSFT %>% stl(s.window = "periodic") %>% seasadj() -> MSFT_SA
autoplot(MSFT_training_SA) + ggtitle("Seasonally Adjusted Microsoft Prices")
#Afit2 <- auto.arima(APPL_training, approximation = TRUE)
# Mfit2 <- auto.arima(MSFT_training, approximation = TRUE)
# summary(Afit2)
# summary(Mfit2)
#
# Afit3 <- auto.arima(APPL_training_SA, approximation = TRUE)
# Mfit3 <- auto.arima(MSFT_training_SA, approximation = TRUE)
# summary(Afit3)
# summary(Mfit3)
#auto.arima not working, so manually use arma() to create model with seasonally adjusted time series
# autoarmafit(APPL_training_SA)
# autoarmafit(MSFT_training_SA)
#ARMA(p,q) models
Afit4 <- arma(APPL_training_SA, order = c(6,3))
summary(Afit4)##
## Call:
## arma(x = APPL_training_SA, order = c(6, 3))
##
## Model:
## ARMA(6,3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.256932 -0.101978 -0.006592 0.089909 6.763310
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.146828 0.182591 -0.804 0.4213
## ar2 -0.504286 0.219250 -2.300 0.0214 *
## ar3 0.435355 0.289374 1.504 0.1325
## ar4 0.038023 0.022003 1.728 0.0840 .
## ar5 -0.031443 0.014042 -2.239 0.0251 *
## ar6 0.063026 0.034537 1.825 0.0680 .
## ma1 0.178053 0.185186 0.961 0.3363
## ma2 0.510486 0.225932 2.259 0.0239 *
## ma3 -0.444735 0.298101 -1.492 0.1357
## intercept 0.013614 0.009157 1.487 0.1371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Fit:
## sigma^2 estimated as 0.223, Conditional Sum-of-Squares = 1598.45, AIC = 9614.88Mfit4 <- arma(MSFT_training_SA, order = c(3,2))
summary(Mfit4)##
## Call:
## arma(x = MSFT_training_SA, order = c(3, 2))
##
## Model:
## ARMA(3,2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.121622 -0.145233 -0.001371 0.145240 4.716149
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.186631 0.030097 6.201 5.61e-10 ***
## ar2 -0.880050 0.046131 -19.077 < 2e-16 ***
## ar3 -0.060310 0.012424 -4.854 1.21e-06 ***
## ma1 -0.224218 0.027813 -8.062 6.66e-16 ***
## ma2 0.906097 0.044303 20.452 < 2e-16 ***
## intercept 0.008690 0.009623 0.903 0.367
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Fit:
## sigma^2 estimated as 0.2352, Conditional Sum-of-Squares = 1686.3, AIC = 9987.83#the seasonally adjusted APPL stocks are best modeled by an ARMA(5,3) with AR(2), AR(3), AR(5), MA(2), and MA(3) being significant
#the seasadj MSFT stocks are best modeled by an ARMA(3,2) with AR(3) and MA(1) being significant
#plot models on top of seasonaly adj series
autoplot(APPL_SA) + ggtitle("Seasonally Adjusted Apple Stock Closing Prices(log-difference)") + autolayer(Afit4$fitted.values, color = "blue")## Warning: Removed 6 row(s) containing missing values (geom_path).
autoplot(MSFT_SA) + ggtitle("Seasonally Adjusted Microsoft Stock Closing Prices(log-difference)") + autolayer(Mfit4$fitted.values, color = "blue")## Warning: Removed 3 row(s) containing missing values (geom_path).
#plot 1500 step ahead forecast for both models
plot(forecast(Afit4$fitted.values, h = 1500), lwd = 0.25, main = "ARMA(5,3) Forecast for APPL Stock(Seasonally Adjusted)")
plot(forecast(Mfit4$fitted.values, h = 1500), lwd = 0.25, main = "ARMA(3,2) Forecast for MSFT Stock(Seasonally Adjusted)")
# plot(forecast(bfit1,h=20),shadecols="oldstyle", col = "blue")
# grid()
# nberShade()
# lines(MSFTts)
# lines(bfit1$fitted, col = "blue")
#
# plot(forecast(efit1,h=20),shadecols="oldstyle", col = "blue")
# grid()
# nberShade()
# lines(MSFT)
# lines(efit1$fitted, col = "blue")
#
# #look at residuals
# tsdisplay(bfit1$residuals)
# tsdisplay(efit1$residuals)#ets model
Afit5 <- ets(APPL_training_SA)## Warning in ets(APPL_training_SA): I can't handle data with frequency greater
## than 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.summary(Afit5) ## ETS(A,N,N)
##
## Call:
## ets(y = APPL_training_SA)
##
## Smoothing parameters:
## alpha = 0.001
##
## Initial states:
## l = 9e-04
##
## sigma: 0.4753
##
## AIC AICc BIC
## 53032.22 53032.22 53052.85
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.006934227 0.4752225 0.2268962 99.26095 102.4022 0.7963111
## ACF1
## Training set 0.01235755Mfit5 <- ets(MSFT_training_SA)## Warning in ets(MSFT_training_SA): I can't handle data with frequency greater
## than 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.summary(Mfit5)## ETS(A,N,N)
##
## Call:
## ets(y = MSFT_training_SA)
##
## Smoothing parameters:
## alpha = 1e-04
##
## Initial states:
## l = 0.0059
##
## sigma: 0.4862
##
## AIC AICc BIC
## 53357.59 53357.59 53378.22
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0005340325 0.4861207 0.2844913 94.83416 110.7957 0.7120463
## ACF1
## Training set -0.0303318plot(Afit5)
accuracy(Afit5)## ME RMSE MAE MPE MAPE MASE
## Training set 0.006934227 0.4752225 0.2268962 99.26095 102.4022 0.7963111
## ACF1
## Training set 0.01235755plot(forecast(Afit5,level=c(50,80,95)), shadecols="oldstyle", col = "blue")
grid()
nberShade()
lines(APPL_training_SA)
lines(fitted(Afit5), col = "blue")
plot(Mfit5)
accuracy(Mfit5)## ME RMSE MAE MPE MAPE MASE
## Training set 0.0005340325 0.4861207 0.2844913 94.83416 110.7957 0.7120463
## ACF1
## Training set -0.0303318plot(forecast(Mfit5,level=c(50,80,95)), shadecols="oldstyle", col = "blue")
grid()
nberShade()
lines(MSFT_training_SA)
lines(fitted(Mfit5), col = "blue")
# #ARIMA with Linear Trend
# bfit3 = Arima(APPL_training_SA,order=c(0,0,1),xreg = cbind(t),seasonal=list(order=c(0,0,1)))
# summary(bfit3)
# plot(forecast(bfit3$fitted,h=20),shadecols="oldstyle", col = "blue")
# grid()
# nberShade()
# lines(fitted(bfit3), ,col="blue")
# lines(bondts)#HW model
#use on non-transformed data
Afit6 <- HoltWinters(appl_ts)
summary(Afit6)## Length Class Mode
## fitted 33500 mts numeric
## x 8627 ts numeric
## alpha 1 -none- numeric
## beta 1 -none- numeric
## gamma 1 -none- numeric
## coefficients 254 -none- numeric
## seasonal 1 -none- character
## SSE 1 -none- numeric
## call 2 -none- callaccuracy(forecast(Afit6))## ME RMSE MAE MPE MAPE MASE
## Training set 0.02044225 1.489885 0.5278037 -0.02143473 2.276895 0.04999679
## ACF1
## Training set -0.01346426Mfit6 <- HoltWinters(msft_ts)
summary(Mfit6)## Length Class Mode
## fitted 33500 mts numeric
## x 8627 ts numeric
## alpha 1 -none- numeric
## beta 1 -none- numeric
## gamma 1 -none- numeric
## coefficients 254 -none- numeric
## seasonal 1 -none- character
## SSE 1 -none- numeric
## call 2 -none- callaccuracy(forecast(Mfit6))## ME RMSE MAE MPE MAPE MASE
## Training set 0.01182961 0.88616 0.406883 0.006770798 1.533883 0.05875956
## ACF1
## Training set -0.02139822#try HW model transformed data
Afit7 <- HoltWinters(APPL_training)
summary(Afit7)## Length Class Mode
## fitted 27692 mts numeric
## x 7175 ts numeric
## alpha 1 -none- numeric
## beta 1 -none- numeric
## gamma 1 -none- numeric
## coefficients 254 -none- numeric
## seasonal 1 -none- character
## SSE 1 -none- numeric
## call 2 -none- callaccuracy(forecast(Afit7))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.004320891 0.4936073 0.2104607 NaN Inf 0.7386293 0.01367678Mfit7 <- HoltWinters(MSFT_training)
summary(Mfit7)## Length Class Mode
## fitted 27692 mts numeric
## x 7175 ts numeric
## alpha 1 -none- numeric
## beta 1 -none- numeric
## gamma 1 -none- numeric
## coefficients 254 -none- numeric
## seasonal 1 -none- character
## SSE 1 -none- numeric
## call 2 -none- callaccuracy(forecast(Mfit7))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01102449 0.5048729 0.2784849 NaN Inf 0.6970131 -0.02868174plot(Afit6)
plot(Mfit6)
plot(forecast(Afit6,level=c(50,80,95)))
grid()
nberShade()
plot(forecast(Mfit6,level=c(50,80,95)))
grid()
nberShade()
plot(Afit7)
plot(Mfit7)
plot(forecast(Afit7,level=c(50,80,95)))
grid()
nberShade()
plot(forecast(Mfit7,level=c(50,80,95)))
grid()
nberShade()
#HW additive seasonal models for APPL
# Afit6 <- hw(APPL_training,seasonal="additive")
# Afit7<- hw(APPL_training,seasonal="multiplicative")
# autoplot(APPL_training) +
# autolayer(Afit6, series="HW additive forecasts", PI=FALSE) +
# autolayer(Afit7, series="HW multiplicative forecasts",
# PI=FALSE) +
# xlab("Year") +
# ylab("Log-Difference of APPL Stock Closing Prices") +
# ggtitle("HW Additive") + guides(colour=guide_legend(title="Forecast"))
# #
# #HW additive seasonal models for MSFT
# Mfit6 <- HoltWinters(MSFT_training,seasonal="additive")
# #Mfit7<- HoltWinters(MSFT_training,seasonal="multiplicative")
# autoplot(MSFT_training) +
# autolayer(Mfit6, series="HW additive forecasts", PI=FALSE) +
# # autolayer(Mfit7, series="HW multiplicative forecasts",
# # PI=FALSE) +
# xlab("Year") +
# ylab("Log-Difference of MSFT Stock Closing Prices") +
# ggtitle("HW Additive") +
# guides(colour=guide_legend(title="Forecast"))Part (c)
We are going to move forward using our seasonally adjusted ARMA models. When we plot the residuals vs the fitted values of our seasonally adjusted ARMA(5,3) model for APPL Stocks and our seasonally adjusted ARMA(3,2) for MSFT Stocks, we see that the residuals look very good. For both models, they are randomly and very densely clustered around zero. In both cases, our range of values for our residuals is very small as well which is also a good sign.
#plot the residuals vs fitted values for our seasadj ARMA models
plot(Afit4$residuals, Afit4$fitted, main = "Seasonally Adjusted ARMA(5,3) Residuals vs Fitted Values")
plot(Mfit4$residuals, Mfit4$fitted, main = "Seasonally Adjusted ARMA(3,2) Residuals vs Fitted Values")
Part (e)
Our residuals from our ARMA models do show significant lag spikes throughout and slowly decay towards 0. However, the ACF & PACF otherwise look very good and show little dynamics are left - a sign of WN. We use our Box tests to see if our residuals are WN and indeed, we get high p-values for all around and can conclude that our residuals are WN and our model has little unexplained dynamics left.
#use tsdisplay to plot the ACF & PACF of the residuals for the model
tsdisplay(Afit4$residuals)
tsdisplay(Mfit4$residuals)
# Test for White noise using Ljung-Box and Box-Pierce
# Ljung-Box Q-Statistic
Box.test(Afit4$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Afit4$residuals
## X-squared = 0.68196, df = 1, p-value = 0.4089# Box-Pierce Q-Statistic
Box.test(Afit4$residuals, type = "Box-Pierce")##
## Box-Pierce test
##
## data: Afit4$residuals
## X-squared = 0.68168, df = 1, p-value = 0.409#low p val regects that data is not WN so need to continue modeling data
# Ljung-Box Q-Statistic
Box.test(Mfit4$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Mfit4$residuals
## X-squared = 0.34981, df = 1, p-value = 0.5542# Box-Pierce Q-Statistic
Box.test(Mfit4$residuals, type = "Box-Pierce")##
## Box-Pierce test
##
## data: Mfit4$residuals
## X-squared = 0.34966, df = 1, p-value = 0.5543Part (f)
#define recursive residuals
Arec=recresid(Afit4$residuals~1)
plot(Arec, pch=16,ylab="Recursive Residuals")
plot(recresid(Afit4$residuals~1, type = "Rec-CUSUM"), main = "Rec-CUSUM Plot Afit4 rec res")## Warning: In lm.fit(x[1:q, , drop = FALSE], y1, tol = qr.tol, ...) :
## extra argument 'type' will be disregarded## Warning: In lm.fit(x[1:(r - 1), , drop = FALSE], y1, tol = qr.tol, ...) :
## extra argument 'type' will be disregarded
Mrec=recresid(Mfit4$residuals~1)
plot(Mrec, pch=16,ylab="Recursive Residuals")
plot(recresid(Mfit4$residuals~1, type = "Rec-CUSUM"), main = "Rec-CUSUM Plot Mfit4 rec res")## Warning: In lm.fit(x[1:q, , drop = FALSE], y1, tol = qr.tol, ...) :
## extra argument 'type' will be disregarded
## Warning: In lm.fit(x[1:(r - 1), , drop = FALSE], y1, tol = qr.tol, ...) :
## extra argument 'type' will be disregarded
Part (i)
Now, we will fit a VAR model to our data and see if we can determine if there is a deterministic relationship between the two variable. Its hard to see
#
plot(appl_ts, main = "Apple Vs. Microsoft Historical Stock Prices")
nberShade()
grid()
lines(msft_ts, col = "Blue")
lines(appl_ts)
plot(appl_ts, main = "Apple Vs. Microsoft Historical Stock Prices", xlim = c(2008,2013), ylim = c(-5, 150))
nberShade()
grid()
lines(msft_ts, col = "Blue")
lines(appl_ts)
plot(APPL, main = "Apple Vs. Microsoft Historical Stock Prices(d=1)", xlim = c(2009.4,2009.75))
nberShade()
grid()
lines(msft_ts, col = "Blue")
lines(APPL)
plot(APPL, main = "Apple Vs. Microsoft Historical Stock Prices(d=1)", xlim = c(2009.47,2009.56), ylim = c(-2, 2))
nberShade()
grid()
lines(MSFT, col = "Blue")
lines(APPL)
ccf(appl_ts,msft_ts,ylab="Cross-Correlation Function", main = "APPL & MSFT STOCK CCF")
ccf(APPL,MSFT,ylab="Cross-Correlation Function", main = "APPL & MSFT STOCK CCF")
#use varselect to see what lag order to use
#we need to combine the variables into 1 data frame first:
# y1=cbind(appl_ts, msft_ts)
# y2=cbind(APPL_training, MSFT_training)
# #y_ts=ts.union(starts, comps) # You can also use this function
# y1_tot=data.frame(y1)
# y2_tot=data.frame(y2)
# head(y1_tot)
# head(y2_tot)
#
# VARselect(y1_tot, lag.max = 252)
# VARselect(y2_tot, lag.max = 252)Citations:
[0] Macrotrends Stock Research Data. (2020). Microsoft - 34 Year stock Price History. Retrieved from https://www.macrotrends.net/stocks/charts/MSFT/microsoft/stock-price-history [1] Macrotrends Stock Research Data. (2020). Apple - 40 Year Stock Price History. Retrieved from https://www.macrotrends.net/stocks/charts/AAPL/apple/stock-price-history