1
$\begingroup$

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:

enter image description here

I checked the ACF and PACF plots and since the data wasn't stationary I transformed the data with diff.

enter image description here

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)[12] 

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)[12]
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).

enter image description here

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.

enter image description here

$\endgroup$
1
$\begingroup$

I am not a statistician so what I say is based on what I have read in various statistical sources, including at this site. These comments are general not related specifically how to improve your model. I would try using auto.arima to estimate the ARIMA/SARIMA. The classical approach you try supposedly does not work well with mixed P and Q (or may not work well anyhow). The rule of thumb is you need 50 monthly data points at a minimum probably more. My guess is 4 months, your hold out sample is also too little to test your model, I use 12 although I know of no rule of thumb. Just because you difference the data once does not make it stationary. I would test the differenced data for stationarity. The various stationarity tests have weak power, so look at the data itself for trend and use multiple tests of stationarity. There are many reasons you may not be able to predict the future. Things can change between the past and the future where you estimate results. ARIMA makes assumptions such as linearity that may not be true. You have to remove all serial correlations and it appears that you have not based on your box ljung test. You might want to consider exponential smoothing in addition to ARIMA and possibly mix various models - which has been shown to improve results. I run 6 ESM models, am going to add an ARIMA model to this, and pick the best three to average as result of a 12 month hold out sample. Even then I don't expect it will always work, there are to many uncertainties.

$\endgroup$
5
  • $\begingroup$ Thank you for your comment. I'm using auto.arima. Regarding the data, I imagined the historical data is too short but I don't have more data available. I only need to forecast 4 or so months. I tested the data with difference equal to 1 and it was okay. Hence using d = 1, but I only used the adf.test. My understanding is that auto.arima also doesn't work for d higher than 2. Did I misunderstood the Ljung box test? No significant is white noise so the model should be fine. No correlated residuals,right? $\endgroup$
    – psoares
    Aug 10 '20 at 21:27
  • $\begingroup$ I can't remember if auto.arima will work with d greater than 2 but you are not going to have many series with a d greater than 2 anyway. My recommendation given power limits on adf is to inspect the data to determine if it has a trend and run kpss which has the opposite null of adf. If they agree its more evidence you have stationarity. I thought you said the Boxljung was significant "since the Ljung-Box is significant." The problem with too little data, among other things, is its harder to capture seasonality I believe. $\endgroup$
    – user54285
    Aug 11 '20 at 0:49
  • $\begingroup$ The exponential smoothing was significant for the Ljung box test, but accuracy and the fitted data looks better with the exponential smoothing. Hence my confusion. ARIMA has these weird peaks in the fitted data (last plot, red line). Also, I'm not sure if d equal to 1 is enough to create stationarity and I thought ets also worked with non-stationary data. Thank you! I will check kpss. $\endgroup$
    – psoares
    Aug 11 '20 at 4:59
  • $\begingroup$ I tried kpss and the data is stationary. Even if I remove d = 1 auto.arima adds the differencing so I guess it's appropriate. I also used your suggestion and did the training with three years and the test with one year and accuracy is much better for arima now. I guess between the two it would be best to use arima. I need to look into combine models. How would you combine the CI? Just average as well, or take the lower and upper? ARIMA gives me wide CIs but it's probably to be expected since the historical data is only 4 years. $\endgroup$
    – psoares
    Aug 11 '20 at 8:34
  • $\begingroup$ I don't use the CI, I am running a practitioner not academic exercise. I just average the predictions. For decades no one realized that exponential smoothing (which operates very differently than ARIMA) even had CI. I am not at all sure how one would average CI. I suspect you can make auto.arima not difference by changing the defaults, but I would not unless you had strong reason to do so. $\endgroup$
    – user54285
    Aug 11 '20 at 18:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.