I have the number of hospitalizations due to a particular disease and I'm trying to forecast the number of hospitalizations for the next 3 months.
I've been reading about ARIMA and time series but I'm mainly learning on my own and I have some questions and I would also like to check what I've done is correct. First I remove the known last 4 values to compare the forecasted values with the observed values.
The plot of observed values:
I checked the ACF and PACF plots and since the data wasn't stationary I transformed the data with
Then I proceeded to the model. I fitted an ARIMA with seasonal pattern and differencing equal to 1.
> fit_sarima <- auto.arima(ts, d = 1, seasonal = T, stepwise = F, approximation = F) > summary(fit_sarima) Series: ts ARIMA(0,1,0)(1,1,0) Coefficients: sar1 -0.5835 s.e. 0.1413 sigma^2 estimated as 31949: log likelihood=-219.95 AIC=443.9 AICc=444.3 BIC=446.9 Training set error measures: ME RMSE MAE MPE MAPE Training set -24.67867 149.0825 97.59303 -2.522906 8.812427 MASE ACF1 Training set 0.5977595 -0.05551552 > checkresiduals(fit_sarima) Ljung-Box test data: Residuals from ARIMA(0,1,0)(1,1,0) Q* = 8.9865, df = 8, p-value = 0.3434 Model df: 1. Total lags used: 9
The model seems to fit the data. I forecasted the next 4 months and compared them with the observed values.
> forecast(fit_sarima, h=4, level=95) Point Forecast Lo 95 Hi 95 Nov 2019 823.9865 473.6554 1174.318 Dec 2019 1092.5124 597.0694 1587.955 Jan 2020 1578.2382 971.4469 2185.029 Feb 2020 1478.9154 778.2532 2179.578 > obs time hosp 1: nov 2019 863 2: dez 2019 1029 3: jan 2020 1186 4: fev 2020 1038
However, when I plot the
fitted values with the observed values I get some weird values for the forecast. Maybe it's expected, I'm not expecting a perfect fit, but I was expecting a somewhat more smoothed curve. Should I have more data available to forecast? (I don't have more but I'm assuming this is a limitation for forecast).
I also used a
ets model to fit the data and although the model seemed less appropriate since the Ljung-Box is significant and the MPE and AIC are higher the curve is smoother compared with the ARIMA, and the confidence intervals smaller.
> fit_exp <- ets(ts) > summary(fit_exp) ETS(M,N,M) Call: ets(y = ts) Smoothing parameters: alpha = 0.5077 gamma = 3e-04 Initial states: l = 1156.097 s = 1.1451 0.8585 0.832 0.7256 0.8352 0.8213 0.848 0.9036 1.0145 1.1833 1.3147 1.5183 sigma: 0.1088 AIC AICc BIC 629.7517 645.7517 657.1813 Training set error measures: ME RMSE MAE MPE MAPE Training set -8.479541 127.0509 84.58677 -1.216796 7.032338 MASE ACF1 Training set 0.5180958 0.1998706 > checkresiduals(fit_exp) Ljung-Box test data: Residuals from ETS(M,N,M) Q* = 27.678, df = 3, p-value = 4.244e-06 Model df: 14. Total lags used: 17 > forecast(fit_exp, h=4, level=95) Point Forecast Lo 95 Hi 95 Nov 2019 874.5813 688.0925 1061.070 Dec 2019 1166.6136 887.2940 1445.933 Jan 2020 1546.8097 1139.8502 1953.769 Feb 2020 1339.2916 957.7705 1720.813 > obs time hosp 1: nov 2019 863 2: dez 2019 1029 3: jan 2020 1186 4: fev 2020 1038
Again, looking at the fitted values the ARIMA (red) looks a bit strange. Is there any way to improve this model? Or should I go with the exponential smoothing model instead (blue)?
**I ended up comparing the accuracy of both models considering the 3 extra months the test dataset and the remaining the training dataset. The accuracy was better for the exponential smoothing and I ended up selecting that model instead, despite the fit of the ARIMA being better.