Interrupted time series analysis of data with deterministic trend and seasonality

Question

I am trying to evaluate the impact of an intervention on a selected outcome variable using interrupted time series data. I have aggregated a five-year data into monthly values to create a data-set 0f 60 time periods.
I have checked for serial correlation, normality, trend, and seasonality using Durbin-Watson test, Augmented Dickey-Fuller test, acf and pacf plots of the preliminary OLS model, ad.test, MannKendall test, kpps test, seasdum tests, etc and found the data to be non-serially correlated and normally distributed.
However, since there are deterministic trend and seasonality, as confirmed by MnnKendall and seasdum tests, I cannot consider it to be stationary. The auto.arima function in R also suggests to take ARIMA(0,1,1) model with first order differencing.
I do not want to make predictions using ARIMA models. I want to fit appropraite regression model to determine the impact of the intervention on the selected outcome. A linear regression model is not appropriate for non-stationary data and corCLASSES of the generalized least square function of R does not include ARIMA models; it includes only corARMA (p,q).
How can I examine the effect of the intervention?

I have setup the data for interrupted time series analysis as follows:

1. outcome variable- is continuous that is to be evaluated for change in level and trend following the intervention
   2. Trend is set up to estimate the change in trend following the intervention
   3. Time- Time (in months) (1-24- pre-intervention, 25-60- post-intervention)
   4. Intervention- dichotomous variable to estimate the change in level on the outcome variable following the intervention (0- pre-intervention; 1- post-intervention)
   5. nur- continuous, independent variable that we thought might contribute to the changes in the outcome variable

Here is the data:

outcome:

22.07110439
18.90095847
18.91359773
19.20793651
20.88607595
18.20779221
19.15797788
21.25701625
20.88034188
21.86503067
24.01552795
19.9045045
18.5278174
18.46753247
20.83765348
17.94267516
18.60979463
19.90101892
16.67349927
17.83104396
18.35260116
17.97080292
17.26436782
18.6896
21.21225071
19.06471495
17.41436464
16.96333333
16.82025678
18.2585034
19.05191257
15.92080537
17.07025411
16.67902996
17.398017
17.45731707
17.47838617
17.12740741
15.23802612
16.64904552
16.50793651
20.46721311
15.08291457
17.80658436
14.71666667
16.48233861
15.48533724
13.4244898
16.81895504
14.59757739
13.5566879
14.08549223
13.14965035
13.53341289
15.83421751
13.20208605
16.10674157
16.70136986
14.79861111
14.95229983
Trend:-
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

Intervention:

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

Time:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

Nur:

0.571666667
0.57
0.571666667
0.586666667
0.591666667
0.591666667
0.565
0.568333333
0.566666667
0.58
0.571666667
0.57
0.563333333
0.561666667
0.561666667
0.613333333
0.615
0.623333333
0.606666667
0.613333333
0.611666667
0.611666667
0.615
0.611666667
0.625
0.611666667
0.623333333
0.706666667
0.72
0.73
0.716666667
0.715
0.72
0.715
0.71
0.703333333
0.703333333
0.701666667
0.693333333
0.781666667
0.79
0.781666667
0.788333333
0.778333333
0.775
0.775
0.776666667
0.781666667
0.781666667
0.773333333
0.781666667
0.851666667
0.846666667
0.853333333
0.838333333
0.84
0.826666667
0.825
0.825
0.815

Answer 1

My answer to you is in two parts. In the first part I set up AUTOBOX to estimate the model that you have pre=specified. The 60x4 DATA MATRIX is partially shown here

My approach is to simply estimate your model and report the results. You are assuming two intervention variables , one a step and the other a deterministic time trend. The problem with this is that your are not forming the intervention variables analytically ... apparently to eye-balled the original plot of the outcome series and divined the two deterministic predictors . Here is the plot of the Y series called "LIKKA-outcome"

I estimated your model and did not perform any tests of significance for stepping down NOR did I examine the residuals for purposes of model augmentation. Although you stated that you were not interested in forecasting , the pattern of the forecasts speaks volumes about the model and it's adequacy. In addition the residual plot tells a lot about the sufficiency of the assumed model.

Here is the Actual/Fit and Forecast for the model you specified.

The plot of the residuals (clearly non-random )

The acf of the residuals showing no evidence of structure ( due to the untreated anomalies ! )

The model ( presented twice )

and

The user-specified two intervention.deterministic series are not-significant while the NUR SERIES ) the user-suggested stochastic series is significant. In summary your two guesses at the form of the latent deterministic structure are of no value.

SECOND ANALYSIS (using AUTOBOX's expert system):

Now I took your data set using the LIKKA-Outcome (Y) series and the user-suggested stochastic input series (MUR) . . I employed the standard paradigm http://www.autobox.com/pdfs/A.pdf to identify/form/revise a Transfer Function Model and obtained the following model

and here with acf here

Note that a seasonal ar(12) was found along with a sophisticated memory model for the MUR series along with 2 downwards level shifts and a seasonal pulse at period 25 and 4 one-time anomalies. Following is the Actual and CLeansed graph

The Actual/Fit and Forecast plot is here . The plot of the residuals is here

The available r stuff doesn't meet your needs BUT there is an R VERSION of the software that I used (AUTOBOX) which does.