Load Libraries

library(quantmod)
library(MASS)
library(tseries)
library(forecast)

Here we’re gonna be using stock price data for Netflix, Inc. (NFLX) from January 2015 all the way up to December 2019.

getSymbols("NFLX", from="2015-01-01", src = 'yahoo', periodicity = 'monthly')
## [1] "NFLX"
stock_ts <- ts(NFLX$NFLX.Close, start=2015, frequency=12)
plot(stock_ts,type = 'l')
title('NFLX Price')

So, we have 5-year monthly data for Netflix and we’re gonna be using the ARIMA model to conduct price forecasts on this data. We’re gonna do two things.

The first thing is we’re going to get our data ready to run the ARIMA model and then secondly we’re actually going to generate forecasts using the ARIMA model and validate our ARIMA to make sure that we’re getting accurate forecasts.

I’m going to be using the first forty eight data points in my time series and I’m gonna be converting it into logarithmic format. The reason using the first forty eight data points is that is gonna be my training data. So that, I’m gonna to actually build the ARIMA model and then the remainder of the data are going to be compare with the forecasts generated by the ARIMA, so we’re gonna be seeing how accurate is the ARIMA model in predicting the actual prices of Netflix. The reason we want to convert into logarithmic format is because stock prices are based on returns, and returns are based on percentages. So, we’re going to be conversing that in logarithmic format to make it suitable for ARIMA model.

lnstock=log(stock_ts[1:48])
lnstock
##  [1] 4.144947 4.217215 4.086432 4.375757 4.490336 4.541683 4.738914
##  [8] 4.745193 4.637250 4.685644 4.814864 4.739526 4.520048 4.536998
## [15] 4.627225 4.500143 4.630545 4.516120 4.513603 4.579339 4.590564
## [22] 4.827273 4.762174 4.818667 4.946701 4.956742 4.995928 5.025195
## [29] 5.094180 5.006694 5.202137 5.163127 5.200429 5.280306 5.234205
## [36] 5.257287 5.599532 5.674628 5.688161 5.744476 5.862494 5.969807
## [43] 5.821417 5.907213 5.924603 5.709698 5.656446 5.589718

Autocorrelation , Partial Correlation and Phillips-Perron Test

The next thing we’re going to do is I’m going to run auto correlation and partial autocorrelation plots. Auto correlation is the correlation between a data point and its previous. So it’s the correlation between the lags essentially.

acf(lnstock, lag.max = 20)

You can see that we have a very gradual descent in our lags so each line going up represents the degree of correlation between our lags.

pacf(lnstock, lag.max = 20)

And when look at the partial autocorrelation function, we can see that we have an immediate drop in the first lag, meaning that all the higher-order autocorrelations are effectively explained by the first lag autocorrelation.

When we work with a time series we want to ensure that we have a stationary model. If we don’t have a stationary model then we have no way of accurately forecasting mean variance and autocorrelation. It’s going to be quite difficult to infer proper forecasts from that data.

I’m going to get the differences between each time series points and I’m going to be using the Phillips-Perron test to screen for stationarity across both the original time series and the differenced time series.

pp.test(lnstock,alternative="stationary")
## 
##  Phillips-Perron Unit Root Test
## 
## data:  lnstock
## Dickey-Fuller Z(alpha) = -9.3512, Truncation lag parameter = 3,
## p-value = 0.5512
## alternative hypothesis: stationary
difflnstock=diff(lnstock,1)
difflnstock
##  [1]  0.072268007 -0.130782761  0.289324645  0.114579340  0.051346185
##  [6]  0.197231491  0.006278917 -0.107942856  0.048393444  0.129220163
## [11] -0.075337494 -0.219478323  0.016950560  0.090226757 -0.127082269
## [16]  0.130402555 -0.114425090 -0.002517409  0.065736403  0.011224669
## [21]  0.236709156 -0.065099285  0.056493450  0.128033700  0.010041082
## [26]  0.039185485  0.029267774  0.068984174 -0.087485368  0.195442598
## [31] -0.039009334  0.037301402  0.079877199 -0.046100666  0.023081622
## [36]  0.342245355  0.075095865  0.013532840  0.056315318  0.118017751
## [41]  0.107312499 -0.148389295  0.085795589  0.017390368 -0.214905072
## [46] -0.053252001 -0.066728737
pp.test(difflnstock,alternative="stationary")
## 
##  Phillips-Perron Unit Root Test
## 
## data:  difflnstock
## Dickey-Fuller Z(alpha) = -52.964, Truncation lag parameter = 3,
## p-value = 0.01
## alternative hypothesis: stationary

You can see that we have a p-value of 0.5546 for our original ln time series, indicates non stationarity. However for the time series that we differenced, we now have a p-value of 0.01 so that indicates at the 5% level we can reject a null hypothesis of non stationarity.

Auto.arima

What we’re gonna be doing is we’re now going to be running our time series using auto.arima(). We’re going to be running it on the original log configuration because what auto.arima() is going to do is it’s going to choose the number of differences for us automatically. What that’s gonna do is instead of manually choosing the number of auto regression, the order of the moving average model, and the number of differences, auto.arima() is just going to do that for us automatically.

library(forecast)

stock_arima <- auto.arima(lnstock, seasonal = F)
stock_arima
## Series: lnstock 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0307
## s.e.  0.0170
## 
## sigma^2 estimated as 0.01391:  log likelihood=34.29
## AIC=-64.58   AICc=-64.31   BIC=-60.88

We can now see that we have an ARIMA configuration of 0,1,0 with drift.

We come up with our final forecast values and you can see here that we have used our exponent function to convert our logarithms back into into raw price so we’re using the exp function in R and we’re converting our price back into into a proper format.

Forecasted values from ARIMA

forecastedvalues_stock = forecast(stock_arima,h=12)
forecastedvalues_stock
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 49       5.620457 5.469336 5.771579 5.389337 5.851578
## 50       5.651197 5.437479 5.864916 5.324343 5.978051
## 51       5.681937 5.420186 5.943687 5.281624 6.082250
## 52       5.712677 5.410433 6.014920 5.250435 6.174918
## 53       5.743417 5.405498 6.081335 5.226615 6.260218
## 54       5.774156 5.403985 6.144328 5.208028 6.340284
## 55       5.804896 5.405065 6.204727 5.193408 6.416384
## 56       5.835636 5.408199 6.263073 5.181928 6.489344
## 57       5.866376 5.413010 6.319741 5.173013 6.559738
## 58       5.897115 5.419227 6.375004 5.166247 6.627984
## 59       5.927855 5.426641 6.429069 5.161314 6.694396
## 60       5.958595 5.435094 6.482096 5.157969 6.759221
plot(forecastedvalues_stock)

forecastedvaluesextracted = as.numeric(forecastedvalues_stock$mean)
finalforecastedvalues=exp(forecastedvaluesextracted)
finalforecastedvalues
##  [1] 276.0156 284.6320 293.5174 302.6802 312.1290 321.8728 331.9207
##  [8] 342.2823 352.9674 363.9860 375.3486 387.0659

Percentage Error

I’m gonna be taking the remainder of the data comparing it to the forecasts that actually generated.

df<-data.frame(stock_ts[49:60],finalforecastedvalues)
col_headings<-c("Actual Price","Forecasted Price")
names(df)<-col_headings
attach(df)
percentage_error=((df$'Actual Price'-df$'Forecasted Price')/(df$'Actual Price'))
percentage_error
##  [1]  0.18699386  0.20516058  0.17680783  0.18313767  0.09074521
##  [6]  0.12372658 -0.02765013 -0.16521642 -0.31891265 -0.26643480
## [11] -0.19287050 -0.19623555
mean(percentage_error)
## [1] -0.01672903

We can see that our mean percentage error is roughly 1.6% so there’s a roughly 1.6 percent deviation on average between the prices that ARIMA is forecasting for us and the and the actual prices.

Ljung-Box : To check if our residuals are random

Next, we use the Ljung-Box test to see if our residuals are random. The reason to do this is because if there is correlation between our residuals then that’s going to cause problems in our time series model. It’s not good to have correlations between your residuals, that could skew the accuracy of the estimates that we’re getting and our model won’t be useful for forecasting..

We’re going to be running the long box test with our null hypothesis being that our residuals are random. The process for picking the lags here can be quite arbitrary. I’m just going with lag 5, 10, and 15.

Box.test(stock_arima$resid, lag = 5, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  stock_arima$resid
## X-squared = 4.632, df = 5, p-value = 0.4624
Box.test(stock_arima$resid, lag = 10, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  stock_arima$resid
## X-squared = 10.57, df = 10, p-value = 0.392
Box.test(stock_arima$resid, lag = 15, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  stock_arima$resid
## X-squared = 11.761, df = 15, p-value = 0.697

We can see that our p-values are all above 0.05. What this means is that we can’t reject the null hypothesis that our residuals are random and the long box test again provides another degree of evidence that the time series is suitable for forecasting.