Forecast iPhone sales for the next 4 quarters.

First, I visualize iphone salest to check if there is any trend and seasonality. From the graph below, there is an increasing trean and peak in September so there is a seasonalility element that we need to address before the forecast. To address that, I take a log to stablize the variance then take a seasonal difference. After doing a some tests, I pick the model that gives the lowest AICC to forecast sales for next 4 quarters. The best model is ARIMA(2,1,0), seasonal(0,1,0).

library(astsa)
library(fpp)
## Loading required package: forecast
## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'default/
## America/Los_Angeles'
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:astsa':
## 
##     gas
## Loading required package: fma
## 
## Attaching package: 'fma'
## The following objects are masked from 'package:astsa':
## 
##     chicken, sales
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
## 
## Attaching package: 'fpp'
## The following object is masked from 'package:astsa':
## 
##     oil
iphone<-read.csv('iphonesales.csv')
iphone_ts<-ts(iphone$sales, start=c(2007, 3), end=c(2017, 1), frequency=4) 
tsdisplay(iphone_ts) #increasing tread and peak on Sep so it is seasonal

tsdisplay(diff(log(iphone_ts),4))

y<-(diff(log(iphone_ts),4))
adf.test(y)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -2.2923, Lag order = 3, p-value = 0.4596
## alternative hypothesis: stationary
kpss.test(y)
## Warning in kpss.test(y): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 1.377, Truncation lag parameter = 1, p-value = 0.01
auto.arima(iphone_ts,lambda=0) #ARIMA(0,1,3)(1,1,0)[4]  #AICc=30.22
## Series: iphone_ts 
## ARIMA(0,1,3)(1,1,0)[4]                    
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ma1     ma2     ma3     sar1
##       -0.8059  0.1974  0.4222  -0.5228
## s.e.   0.1779  0.2610  0.1753   0.1809
## 
## sigma^2 estimated as 0.1073:  log likelihood=-9.04
## AIC=28.07   AICc=30.22   BIC=35.7
fit1<- Arima(iphone_ts, order=c(0,1,1), seasonal=c(0,1,0),lambda=0)#AICc=35.87 
fit2<- Arima(iphone_ts, order=c(0,1,3), seasonal=c(1,1,0),lambda=0)#AICc=30.22
fit3<- Arima(iphone_ts, order=c(2,1,3), seasonal=c(1,1,0),lambda=0)#AICc=25.5
fit4<- Arima(iphone_ts, order=c(2,1,0), seasonal=c(0,1,0),lambda=0)#AICc=23.47

#Pick the best model with lowest AICC
fit<-  Arima(iphone_ts, order=c(2,1,0), seasonal=c(0,1,0),lambda=0)#AICc=23.47
summary(fit)
## Series: iphone_ts 
## ARIMA(2,1,0)(0,1,0)[4]                    
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ar2
##       -0.9476  -0.7233
## s.e.   0.1187   0.1207
## 
## sigma^2 estimated as 0.09623:  log likelihood=-8.34
## AIC=22.67   AICc=23.47   BIC=27.25
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.666146 9.143921 6.018469 -12.58554 23.87181 0.8450193
##                   ACF1
## Training set 0.5556705
Acf(residuals(fit)) #no autocorrelation in the residual

Box.test(residuals(fit), lag=8, fitdf=2, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(fit)
## X-squared = 5.6849, df = 6, p-value = 0.4594
fcast<-forecast(fit,h=4)
plot(fcast)