Lab6

DSO 522: Applied Time Series for Forecasting Box-Jenkins aka ARMA Models AAA Washington Case In 1993, AAA Washington was one of the two regional automobile clubs affiliated with the American Automobile Association (AAA or Triple A) operating in Washington State. Club research had consistently shown that the emergency road service benefit was the primary reason that people join AAA. Providing emergency road service was also the club’s single largest operating expense. It was projected that delivering emergency road service would cost $9.5 million, 37% of the club’s annual operating budget, in the next fiscal year. Michael DeCoria objective is to find a way to predict emergency road service call volume for future years.The data on emergency road service call volume is given in AAA.csv . You have recommended Michael to model these patterns using ARIMA models. He offered you a contract as a consultant to help him apply Box-Jenkins models. Task: Study the calls volume auto-correlations and build a Box-Jenkins model for the AAA call volume. Your Action Plan: (1.) Visualize your data; (2.) Find the ACF and PACF; Propose several models based on your analysis of ACF and PACF. (3.) Choose the model that is the most accurate, robust, and the simplest. (4.) Test the assumptions (5.) Provide recommendations to Mr. DeCoria.

(1.) Visualize your data Use the monthly data of the calls volume, Calls , is recorded from May, 1988 till December, 1992. Produce a time plot of the data, label the axis nicely. Comment on the patterns present in the time series. Include your script. Comment on the scatterplot. data = read.csv ( 'AAA.csv' , header = T) data.ts = ts (data[, 3 ], start = c ( 1988 , 5 ), frequency = 12 ) plot (data.ts, main= 'AAA monthly call volume' , ylab= "Call volume" ) The plot shows level, noise, and seasonality. Thus, we have two choices. First, de- seasonalize first and fit the model. Second, don’t de-season first. But look into ACF and PACF to see how to pick the number.Need to pay attention here, when we fit the model in the first senerio, don’t process based on the de-season data but still using the original one to do the fair comparison with the latter one. Which means the only difference will be the number and the order in the (). (2.) Split the time series into training and validation sets. Leave the last 12 months for the testing set. Use the training set to get ACF and PACF plots. Based on these charts, comment on the patterns present in the calls volume. Suggest at least 2 Box-Jenkings models. Provide

#The training set shows seasonal part tailed off wave pattern in ACF. PACF shows the first lag cut off (after one long bar no more behind) --> SAR(1). *(1,0,0) #Then we look at the non-seasonal part, which is between the seasonal repeating bar.PACF cut off (drops after the first bar), in ACF can be tail off --> AR(1) (1,0,0)* #or cut off --> ARMA(1,1) or AR(1) + MA(1) dtrain = diff (train.ts, 12 ) par ( mfrow= c ( 2 , 2 )) Acf (dtrain, 36 , main = "" ) Pacf (dtrain, 36 , main = "" ) par ( mfrow= c ( 1 , 1 ))

#no season no trend #no-seasonal part, not sure which one better. Can try different combination (1,0,1)*(0,1,0); (1,0,0)*(0,1,0); (0,0,1)*(0,1,0) (3.) Fit the suggested 2 models and access their accuracy, robustness, assumptions. Visual analysis techniques suffice. Include your script, comments. Report all supporting material. m1 = arima (train.ts, order = c ( 0 , 0 , 0 ), seasonal = list ( order= c ( 0 , 1 , 0 ), period = 12 )) summary (m1) ## ## Call: ## arima(x = train.ts, order = c(0, 0, 0), seasonal = list(order = c(0, 1, 0), ## period = 12)) ## ## ## sigma^2 estimated as 2251271: log likelihood = -279.44, aic = 560.88 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE

m2 = arima (train.ts, order = c ( 1 , 0 , 0 ), seasonal = list ( order= c ( 1 , 0 , 0 ), period = 12 )) summary (m2) ## ## Call: ## arima(x = train.ts, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), ## period = 12)) ## ## Coefficients: ## ar1 sar1 intercept ## 0.5002 0.5238 21360.1596 ## s.e. 0.1302 0.1371 582.9503 ## ## sigma^2 estimated as 1390592: log likelihood = -375.7, aic = 759.39 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set -66.93572 1179.233 931.6528 -0.6402688 4.43517 0.7961737 ## ACF1 ## Training set 0.07932037 m2.p = forecast (m2, h= 12 ) accuracy (m2.p, valid.ts) ## ME RMSE MAE MPE MAPE MASE ## Training set -66.93572 1179.233 931.6528 -0.6402688 4.435170 0.8231512 ## Test set 436.35692 1326.482 815.2446 1.6439259 3.572109 0.7203001 ## ACF1 Theil's U ## Training set 0.07932037 NA ## Test set 0.30183881 1.00346 checkresiduals (m2)

## ## Ljung-Box test ## ## data: Residuals from ARIMA(1,0,0)(1,0,0)[12] with non-zero mean ## Q* = 4.6145, df = 7, p-value = 0.7069 ## ## Model df: 2. Total lags used: 9 Compared model 1 and model 2, we can see model 2 has a lower MAPE on testing set, although RMSE is slightly higher. When we look into the residuals, model 2’s residuals are more normally distributed, and no significance shows in ACF. The overall spread is pretty random. In a word, assumptions are met better than model 1. Thus, model 2 is a better model. Propose the best model for the calls volume forecasting. Report its estimated accuracy, comment on its robustness. (4.) Use the best model to produce the forecast for 1993. Visualize your results. m = arima (data.ts, order = c ( 1 , 0 , 0 ), seasonal = list ( order= c ( 1 , 0 , 0 ), period = 12 )) m.pred = forecast (m, h= 12 )

The best model to predict would be SARIMA(1,0,0)*(1,0,0)12. The prediction error should be expected around 3.57%. Based on this model, the average January and February 1993 prediction would be 24055.18 and 22448.91. To better incorporate the budget planning, we should also consider the range between 21764.97 to 26345.39 for January, and between 19854.29 to 25043.52 for February.