# Trying to Refine SARIMA models

I have a SARIMA forecast from statewide Real Estate Sales data.. but I'm not happy with it. The SARIMA parameters are confusing to say the least.

I am finding that the current model is not forecasting high enough, although the month by month fluctuations look reasonable. Currently the projected year ahead is 4.9% above this year. but the year over year gain for the current year is about 10%+, so in short the model is not weighted heavily enough to the recent year over year change. Between the ARIMA and Seasonal components I can't see how to fix/specify this? The Parameter descriptions remain a mystery no matter how many versions I read. Looking at the chart it is clear that the growth rate over the last three years is higher than the current forecast, Surely there is a simple way to deal with this sort of specification in SARIMA, Bonus if it is understandable by mere mortals!?

R code follows (source data is public) I'm a real newb at R as well.

require(RCurl)
fullmatrix <- read.csv(textConnection(myCsv), stringsAsFactors = FALSE)
fullmatrix$date <- as.Date(fullmatrix$date, format = "%m/%d/%y")
fullmatrix$Unit.Sales <- as.numeric(gsub(",", "", fullmatrix$Unit.Sales))
fullmatrix$Average.Selling.Price <- as.numeric(gsub(",", "", fullmatrix$Average.Selling.Price))
fullmatrix$Median.Selling.Price <- as.numeric(gsub(",", "", fullmatrix$Median.Selling.Price))
fullmatrix$Total.For.Sale <- as.numeric(gsub(",", "", fullmatrix$Total.For.Sale))
order.date <- order(fullmatrix$date ) fullmatrix <-fullmatrix[order.date, ] library(xts) require(forecast) require(astsa) sarima(fullmatrix$Unit.Sales,12, 0, 1, 1, 0,1,0) # (slightly higher fore)
fore <- sarima.for(fullmatrix$Unit.Sales,12, 12, 0, 1, 1,0,1,12)#specifies same model fullmatrix.xts <- as.xts(x=fullmatrix[,-1],order.by=fullmatrix$date)
unitsales.xts <- as.xts(as.numeric(fullmatrix.xts$Unit.Sales),order.by=index(fullmatrix.xts)) forecast.xts <- as.xts(as.numeric(fore$pred), order.by=as.Date(c("2013-11-01","2013-12-01","2014-01-01","2014-02-01","2014-03-01","2014-04-01","2014-05-01","2014-06-01","2014-07-01","2014-08-01","2014-09-01","2014-10-01")))
orig_fore <-rbind(unitsales.xts,forecast.xts)
colnames(orig_fore) <- "Unit.Sales"
plot(orig_fore)


• I'm not sure I understand you question. Do you mean that you expected the forecasts to be a little higher, that is, that they follow a slightly increasing pattern with respect to the last years in the sample? Do you expected forecasts with the seasonal pattern that is already forecasted but shifted slightly upwards? I have checked some alternative models for this data, but I do not know if they really address your question. Please clarify a little bit your question. – javlacalle Jun 26 '14 at 17:49
• @javlacalle you are exactly right on my frustration, there does not seem to be a good way that I can see to weight the more recent results to cause the projection to continue to shift upwards.. the seasonal pattern seems to overwhelm the year to year drift upwards.. – dartdog Jun 27 '14 at 3:41

This is what I see when I fit your model. Once the data are loaded

require(RCurl)
require(astsa)
m <- read.csv(textConnection(myCsv), stringsAsFactors = FALSE)
x <- as.numeric(gsub(",", "", as.vector(m[,"Unit.Sales"])))
x <- ts(rev(x), start = c(2000, 5), frequency = 12)
x0sub <- window(x, start=c(2013, 11))
x <- window(x, end = c(2013, 10))


fit <- sarima(x, 12, 0, 1, 1, 0, 1, 12)
fore <- sarima.for(x, 12, 12, 0, 1, 1, 0, 1, 12)
fit
# [trimmed]
# Coefficients:
#          ar1      ar2     ar3      ar4     ar5      ar6     ar7      ar8
#       1.7031  -0.6764  0.1431  -0.4395  0.5119  -0.4085  0.2630  -0.2326
# s.e.  0.0766   0.1492  0.1532   0.1483  0.1534   0.1593  0.1583   0.1542
#          ar9     ar10    ar11     ar12      ma1    sar1     sma1      xmean
#       0.4650  -0.5451  0.3632  -0.1543  -1.0000  0.9999  -0.9797  3773.7099
# s.e.  0.1517   0.1539  0.1525   0.0797   0.0045  0.0001   0.0096    67.2395
# [trimmed]
sum(coef(fit$fit)[1:12]) #[1] 0.9930268 as.vector(fore$pred)
# [1] 3166.736 3309.882 2827.121 3237.666 3989.075 4183.633 4580.797 4665.480
# [9] 4564.829 4529.723 3997.178 3846.385


The sum of the AR coefficients is close to one, it suggests a unit root in the AR polynomial that should be handled taking first differences to the data. Similarly, the seasonal AR coefficient, sar1 = 0.9999 suggests using the seasonal differencing filter.

The ARIMA(0,1,1)(1,0,1) model yields:

fit2 <- sarima(x, 0, 1, 1, 1, 0, 1, 12)
fore2 <- sarima.for(x, 12, 0, 1, 1, 1, 0, 1, 12)
fit2
# [trimmed]
# Coefficients:
#           ma1    sar1     sma1  constant
#       -0.2893  0.9995  -0.9521    5.0690
# s.e.   0.0811  0.0019   0.0970   90.6544
# [trimmed]
as.vector(fore2$pred) # [1] 3188.940 3286.837 2693.416 3088.757 3888.024 3945.766 4397.282 4472.638 # [9] 4293.293 4307.229 3722.071 3554.581  For this model, the optimization algorithm did not converge after 100 iterations and the seasonal AR coefficient still suggests the usage of a seasonal difference. Despite these pitfalls, the forecasts are not that bad (see graphic below) Then I fit an apparently more appropriate model, ARIMA(2,1,1)(2,1,0): fit3 <- sarima(x, 1, 1, 1, 2, 1, 0, 12) fore3 <- sarima.for(x, 12, 1, 1, 1, 2, 1, 0, 12) fit3 # [trimmed] # Coefficients: # ar1 ma1 sar1 sar2 # 0.1728 -0.5422 -0.5505 -0.2367 # s.e. 0.2201 0.1897 0.0840 0.0813 # [trimmed] as.vector(fore3$pred)
# [1] 3313.329 3343.138 2890.382 3233.489 3903.440 4078.577 4461.776 4470.415
# [9] 4459.707 4495.670 3909.286 3796.468


Let's see a plot of one-year-ahead forecasts (sorry for the confusing indentation in the code below).

par(mfrow = c(3, 1), mar = c(2.5, 2, 2, 2))
plot(cbind(x, fore$pred), plot.type = "single", type = "n", xlim = c(2005.417, 2014.750), ylim = c(1890, 6000), main = "ARIMA(12,0,1)(1,0,1)") lines(x, type = "p") lines(x) lines(fore$pred, col = "blue", type = "p")
lines(fore$pred, col = "blue") lines(fore$pred + 1.96 * fore$se, lty = 2, col = "red") lines(fore$pred - 1.96 * fore$se, lty = 2, col = "red") plot(cbind(x, fore2$pred), plot.type = "single", type = "n",
xlim = c(2005.417, 2014.750), ylim = c(1890, 6000),
main = "ARIMA(0,1,1)(1,0,1)")
lines(x, type = "p")
lines(x)
lines(fore2$pred, col = "blue", type = "p") lines(fore2$pred, col = "blue")
lines(fore2$pred + 1.96 * fore2$se, lty = 2, col = "red")
lines(fore2$pred - 1.96 * fore2$se, lty = 2, col = "red")
plot(cbind(x, fore3$pred), plot.type = "single", type = "n", xlim = c(2005.417, 2014.750), ylim = c(1890, 6000), main = "ARIMA(2,1,1)(2,1,0)") lines(x, type = "p") lines(x) lines(fore3$pred, col = "blue", type = "p")
lines(fore3$pred, col = "blue") lines(fore3$pred + 1.96 * fore3$se, lty = 2, col = "red") lines(fore3$pred - 1.96 * fore3$se, lty = 2, col = "red")  The forecasts are similar. The confidence intervals for the models with non-zero order of integration (second and third models) are wider, as it uses to be the case compared to the model with no differencing filter. Now that we have 6 more observations we can use these new data to compare the forecasts from each model. par(mfrow = c(3, 1), mar = c(2.5, 2, 2, 2)) plot(fore$pred, type = "b", col = "red", ylim = range(c(x0sub, fore$pred)), main = "ARIMA(12,0,1)(1,0,1)") lines(x0sub, type = "b", col = "blue") legend("topleft", legend = c("forecasts", "observed values"), col = c("red", "blue"), lty = c(1, 1), bty = "n") plot(fore2$pred, type = "b", col = "red", ylim = range(c(x0sub, fore$pred)), main = "ARIMA(0,1,1)(1,0,1)") lines(x0sub, type = "b", col = "blue") legend("topleft", legend = c("forecasts", "observed values"), col = c("red", "blue"), lty = c(1, 1), bty = "n") plot(fore3$pred, type = "b", col = "red", ylim = range(c(x0sub, fore\$pred)),
main = "ARIMA(2,1,1)(2,1,0)")
lines(x0sub, type = "b", col = "blue")
legend("topleft", legend = c("forecasts", "observed values"),
col = c("red", "blue"), lty = c(1, 1), bty = "n")


The forecasts are actually larger than the new observations in the series and follow the seasonal pattern fairly well. Despite the bad signals that we show in the fit of the second model the forecasts turned to be close to the observed values (this is puzzling to me).

Despite the parameter estimates in some models (e.g. seasonal AR coefficient close to unity) suggested taking differences to the data, it may be a consequence of the level shift around the middle of the sample and the two different trends in the sample (a linear trend in the first sample and more flat in the second part). I tried some other tentative models to capture these facts, the in-sample observations were fitted slightly better by the model but there were no major changes in the forecasts. Just to illustrate this have a look at the autocorrelation functions of the differenced series:

par(mfrow = c(3, 1), mar = c(3.2, 2.5, 3, 2))
acf(diff(x), lag.max = 84, main = "")
mtext("ACF original series", side = 3, adj = 0)
acf(diff(window(x, start = time(x)[89])), lag.max = 84, main = "")
mtext("ACF first subsample (until 2007:09)", side = 3, adj = 0)
acf(diff(window(x, end = time(x)[88])), lag.max = 84, main = "")
mtext("ACF second subsample (from 2007:09)", side = 3, adj=0)
mtext("lag order", side=1, line = 2)


When the whole sample is taken, the autocorrelations of seasonal order are high even for high lag orders, suggesting the use of the seasonal differenced filter. The interpretation of the ACF for the data observed in the two subsamples is the opposite, as they decrease fast and fall within the confidence bands after a few years.