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