I would like to forecast the non-stationary time series, involving several crucial a-priori assumptions following from studying of instances of such series.
I've constructed time-averaged one-point probability distribution function approximated by normal distribution. $$\hat p(x) = \frac{1}{\sqrt{2\pi \sigma^2_{\infty}}} \exp\left(-\frac{x^2}{2\sigma^2_{\infty}}\right)$$ From this point of view, I want the forecast $z_t(l)$ not to exceed this when $l \to \infty$. To put it in other words, variance of $z_t(l)$ must be bounded.
The average two-point probability distribution function $\hat p(x_i,i;x_j,j)$ also has been constructed, which led to identification of autocorrelation function. $\rho(j) \approx A j^{-\alpha} $ provided $0<\alpha<0.5$.
At first, Box-Jenkins identification process led me to $ARIMA(0,1,3)$ model, however
I can't have bounded variance until $d \ne 0$ (which follows from equations for BJ weights $\psi_j$). At the same time, I can't use $d=0$ since initial autocorrelation decreases slowly (which is probably evidence of non-stationarity according BJ). This is the main obstacle to me.
Visually, simulation of $ARIMA(0,1,3)$ does not coincide with behaviour of my samples. And correlations of first difference of the series are in the bad agreement with correlations following from the model.
The analysis of residuals shows significant correlations starting lag 3. This is why my initial statement about $ARIMA(0,1,3)$ is incorrect.
Trying to fit different $ARIMA(p,0,0)$ models, I see that there is significant residual correlations close to the lag $p$ for every $p$. It may assume that I need $ARIMA(\infty,0,q)$ model (as limiting choice), for instance fractional ARIMA.
From [1] I've learned about Fractional $ARIMA(p,d,q)$ models which are $ARIMA(\infty,0,q)$ in effect.
I've not found any GNU R packages with support of missing values for this. Missing values seems to be a kind of challenge.
The publications on fractional ARIMA are quite rare. Are such fractional models really used? Maybe there is a good replacement of ARIMA models for my needs? The forecasting is not my major, I have only pragmatic interest.
From different literature (for instance [2]), I learned that it is practically impossible to decide between fractional ARIMA and models with "level shift". However, I have not found the package for GNU R to fit 'level shift' models.
[1]: Granger, Joyeux.: J. of time series anal. vol. 1 no. 1 1980, p.15
[2]: Grassi, de Magistris.: "When long memory meets the Kalman filter: A comparative study", Computational Statistics and Data Analysis, 2012, in press.
Update: to render my own progress and to answer @IrishStat
My statement about two-point probability distribution is incorrect in general. Constructed in this way function will depend on full series length. So, there is a little to extract from this. At least, parameter named $\alpha$ will depend on full series length.
Lists 2 and 3 also have been updated.
My data is available as dat file here.
At the current moment, I doubt between FARIMA and level shifts, and I still can't find appropriate software to check this options. This is also my first experience with model identification, so any help will be appreciated.
. A significant change point was detected at period 137 suggesting time-varying parameters. The remaining 668 observations suggest a pdq ARIMA Model (3,0,0) with a level.step shift supporting your preliminary conclusions about lag 3. 
. The Actual/Fit/Forecast graph is 
The Residual Plot 
and the acf of the residuals is 
. Since the acf of the residuals shows strong structure at periods 5 and 10 , 
you might further investigate seasonal structure at lag 5. I hope this helps.
