Why does my ARMA forecast get smaller over time?

Question

I am a beginner to time series modeling but I am trying to build an ARMA model to describe a set of 24 observations.

-6.73217
-2.86888
-5.37159
-6.51274
-4.17106
-5.7385
-3.34395
-1.94488
-5.46411
-3.18918
-3.6847
-2.1683
-2.32984
-0.76198
-2.18903
1.094239
-4.8123
0.784199
-1.56708
0.143964
1.131119
2.899746
-0.49872
3.121624

The ACF pattern seems to indicate that this is an AR process.

The PACF shows a significant corrleation at lag 2.

Based on this, I fit a ARMA (2,0,0) model to the data

fitdata <- Arima(tsdata,c(2,0,0))

This resulted in the following model

Series: tsdata ARIMA(2,0,0) with non-zero mean

Coefficients: ar1 ar2 mean 0.1023 0.6691 -2.0834 s.e. 0.1489 0.1652 1.3511

sigma^2 estimated as 3.987: log likelihood=-49.69 AIC=107.39
AICc=109.49 BIC=112.1

I then tried to predict 6 time periods ahead but as you can see the predictions get consistently smaller over time and do not seem to capture the overall trend of the data. Can someone tell me what I'm doing wrong?

pred<-predict(fitdata,n.ahead=6)

EDIT: I should have mentioned earlier that according to a Dickey Fuller test the current time series is stationary.

Augmented Dickey-Fuller Test

data: tsdata Dickey-Fuller = -4.0131, Lag order = 1, p-value = 0.0228 alternative hypothesis: stationary

Answer 1

Any stationary AR(2) process will converge towards the mean. You can see this by yourself by taking any initial values $x_1$ and $x_2$, then developing the future expectations by iteratively plugging in the fitted equation

$$ x_t=0.16x_{t-1}+0.7x_{t-2} $$

(where we set the intercept to $0$). For instance:

set.seed(1)
xx <- rnorm(2)
for ( ii in 3:10 ) xx <- c(xx,.16*xx[ii-1]+0.7*xx[ii-2])
plot(xx,type="l")

More information here.

Incidentally, I would recommend against "rolling your own" ARIMA model. Much better to rely on trusted software, like forecast::auto.arima(). Which indeed finds a drift (i.e., trend) term here, along with a first order integration:

library(forecast)
plot(forecast(auto.arima(tsdata),h=6))

I find this forecast somewhat more convincing than the AR(2) one in the question.

Answer 2

ARIMA modeling procedures that ignore the possibility of anomalies often are quite deficient as suggested by @AdamO here Interrupted Time Series Analysis - ARIMAX for High Frequency Biological Data? ... "The correlogram should be calculated from residuals using a model that controls for intervention administration, otherwise the intervention effects are taken to be Gaussian noise, underestimating the actual autoregressive effect." If one uses the original acf/pacf rather than a acf/pacf conditional upon detected interventions then there often can be poorer results.

I took your 20 values and analyzed them with a robust approach which identified 4 pulses and an ar(1) arima component with only lag 2 .

Note Well that this is a totally different model from what has been represented here suggesting memory of order 2 and order 2 only not .16 and .7 respectfully but o. and .805 and the conditional mean is not 0.0 but -.05268 .

Here is the identified model in two views . and with model statistics here

Some may bridle about using 6 coefficients for 20 observations but they should note that all coefficients are statistically significant yielding a necessary and sufficient model.

The Actual/Fit and Forecast are here with forecasts converging to an assymptote

EDITED after @whubers very insightful and thought-provoking questions.

When I ran it the first time I took the option of simply providing presumptive confidence limits based upon normality of the residuals. Here is the distribution of the residuals (hardly normal).

This time I unleashed the monte-carlo bootstrapping option (2000 realizations per forecast period) where the residuals get re-sampled and limits are then based upon them. Furthermore I optionally enabled pulses to be present in the simulations thus 33.33% ( 4 out of 12 ) where you had correctly guessed 1 of 4 ) had pulses in the realizations because 4 pulses were found in the 12 (18-12) non-zero residuals.. 20 observations minus 2 for the ar and minus 6 for the estimated parameters .In this way a more realistic assumption i.e. the possible re-occurrence of pulses in the future are put in place.

In this manner probability density functions can be constructed based upon the actual distribution of residuals without having to give that distribution a name.

and the forecasts here

Prof. Sam Savage https://www.probabilitymanagement.org/ was inspirational in motivating this particular developmemt of AUTOBOX to meet his simulation/analytical requirements as he wished to have realistic forecasting limits.

A table comparing standard normality-assuming limits(right-hand side of table ) versus the monte-carlo/simulation approach (left) is illuminating.

As usual your thoughtful reflections motivate (TEASE OUT) more critical detail from me.

EDITED AFTER OP'S FIRST COMMENT:

To discuss the differences between the two models , one simply has to compare the model, fit and forecast. I used your ar(2) model form and estimated parameters and obtained a somewhat surprisingly different and better solution in terms of estimated variance . and here with forecasts that are getting smaller as you suggested BUT are better characterized as converging to an assymptote that is lower than recent values and are not truly going to get get smaller forever as intimated by your question.

Your models forecasts are here with Autobox's forecasts here from it's automatically identified model