# How to make this data stationary

What is the correct way to make this data stationary (without cutting it in half)?

Quarterly Data:

17996666000
17395339056
17338377000
17311651000
18043497534
18716063000
19335859000
19300627000
19602120000
19537854000
19888609000
19760257000
19339309125
18885771000
20311922000
19708371000
19326807000
18733706000
21059811000
20523668000
21906648000
24469040000
29176362000
28455903000
30348154000
31588708000
33588602000
32669267000
34233556000
35470378000
36709046000
35698357000
36881680000
37560603000
38169265000
37147402000
38707144000
39120015000
40466332000
41033858000
43142309000



Plot of data :

Differences won't make it stationary

adf.test(data)

Augmented Dickey-Fuller Test

data:  data
Dickey-Fuller = -1.602, Lag order = 3, p-value = 0.7301
alternative hypothesis: stationary

Augmented Dickey-Fuller Test

data:  diff(data)
Dickey-Fuller = -1.4163, Lag order = 3, p-value = 0.8026
alternative hypothesis: stationary

Augmented Dickey-Fuller Test

data:  diff(diff(data))
Dickey-Fuller = -2.6651, Lag order = 3, p-value = 0.3129
alternative hypothesis: stationary

Augmented Dickey-Fuller Test

data:  diff(diff(diff(data)), lag = 4)
Dickey-Fuller = -2.9252, Lag order = 3, p-value = 0.2138
alternative hypothesis: stationary


plot of differenced data:

detrend with a line:

Residuals definitely not stationary:

2 diffs and seasonal difference, still not unit-root stationary per ADF

adf.test(diff(diff(diff(lm(y~t, g)$residuals)), lag=4, differences = 1)) Augmented Dickey-Fuller Test data: diff(diff(diff(lm(y ~ t, g)$residuals)), lag = 4, differences = 1)
Dickey-Fuller = -2.9252, Lag order = 3, p-value = 0.2138
alternative hypothesis: stationary


This data gets stationary if I cut it in half, which eliminates the beginning flat part of the curve. For learning I am interested in seeing how to make it stationary with this part included though, but just for reference:

adf.test(diff(data[17:length(data)]))

Augmented Dickey-Fuller Test

data:  diff(data[17:length(data)])
Dickey-Fuller = -3.8132, Lag order = 2, p-value = 0.03477
alternative hypothesis: stationary


I also tried tsoutliers package to look for interventions, but it gives me errors :

tsoutliers::tso(data,types = c("AO", "TC","LS", "IO"),maxit.iloop=10, maxit.oloop=10)

Error in arima(y, order = fit$$arma[c(1, 6, 2)], seasonal = list(order = fit$$arma[c(3,  :
non-stationary seasonal AR part from CSS

tsoutliers::tso(diff(data),types = c("AO", "TC","LS", "IO"),maxit.iloop=10, maxit.oloop=10)

Error in auto.arima(x = c(-601326943.954987, -56962056.0450134, -26726000,  :
No suitable ARIMA model found
In sqrt(diag(fit\$var.coef)[id]) : NaNs produced


How can I make this stationary without cutting the data in half?

UPDATE:

I performed an intervention analysis following steps from the paper by Tsay, using the model ARIMA(0, 1, 0)(1, 0, 0)4 as an estimate, and following the iterative procedure in the paper, I found the following :

           w            v     hyp   types  t
11 102080630 9.828064e+14 3.256187    IO 22
4   69408316 5.806115e+14 2.880505    TC 22
2   39451617 3.232890e+14 2.194165    AO 22
3   70523673 6.465779e+14 2.773476    LC 22


The only hypothesis that passed the minimum criteria mentioned by Tsay is the innovative outlier at t=23 (23 because the data is differenced, so I add 1)

This is what the transformed data looks like:

The transformed data is mostly stationary after 2 differences:

Also, auto.arima believes this to be the model of the transformed data: ARIMA(0,2,2)(1,0,0)[4]

Differences do make it stationary when a level/step shift indicator is introduced to reflect a deterministic effect at period 23 . No need to segment the data ...just use a comprehensive hybrid model found via http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html .

Here is a useful model using AUTOBOX (a time series forecasting package that I have helped to develop) (0,1,0)(1,0,0)4 with a level shift ( i.e. intercept change at period 23 and a pulse at period 22 .The Actual/Fit and Forecast is here . The Residuals are here with acf here suggesting model sufficiency .

The model includes simple differencing and an intercept adjustment at period 23 along with a simple pulse and a seasonal ar component. Quite simple but powerful and elegant.

Your comment "I also tried tsoutliers package to look for interventions, but it gives me errors " suggests to me possible serious software limitations or a possible "user error" which I can't immediately detect .

The Actual/Cleansed graph is illuminating. while the forecast graph illustrates the 95% prediction limits using monte-carlo resampling.

EDITED TO SHOW EXPLICITELY HOW A FORECAST IS MADE:

The model can be expanded via algebra to present how a prediction can be seen as a "regression type model" . Here is a one-period out forecast computation where differences and autoregressive lag structure is incorporated..

• Thanks again! This is great! So, am I correct that the clean data was produced by detrending with a level change and pulse as regressors? Commented Jan 30, 2020 at 3:51
• Yes ! BUT the phrase "detrending with " should be "incorporating" as the term detrending is appropriate when a predictor series is a trending series. In this case the predictor series is simply a pattern of 22 consecutive 0's and 19 consecutive 1's which does not have a trend. Commented Jan 30, 2020 at 6:27
• I should add ... a pulse is a difference of a step ... a step is a difference of a time trend . When one has a hybrid model that includes a differencing operator ( as in this case) the de facto effect of a level shift is a time trend thus the upwards trending expectation. Commented Jan 30, 2020 at 12:19
• I added a visual of how a forecast can be generated using expanded polynomial structure for this case. It like all forecasts is the weighted sum of both history and any needed deterministic structure. Commented Jan 30, 2020 at 13:10
• You are very welcome . Commented Jan 31, 2020 at 0:41

Well, when you're detrending the data, obviously the linear trend is not sufficient. The trend basically looks like a flat line followed by linear increase starting around $$t=20$$; that is not well approximately by a linear trend because it overestimates near the break point. You might want to use a more general parameterization for the trend- loess could go the job. Something to think about.

• I tried loess and it had the same issue. The loess residuals looked quadratic. I didn't include because post was long enough already. Commented Jan 30, 2020 at 3:40