library(TSA)
data(airmiles)
?airmiles
plot((airmiles),ylab='(airmiles)',xlab='Year', )

plot(log(airmiles),ylab='Log(airmiles)',xlab='Year', )

?filter
S=c(rep(0,2),rep(1,108))
# step pulse of the post-intervention period
plot(ts(S),type='p')
abline(h=0)

omega=-1
delta=.9
m=filter(x=omega*S, filter=c(delta),method='recursive')
plot(ts(m),type='p',main='change in the mean function due to intervention at T=3')

#identify the model for the natural process over the pre-intervention period
plot(window(log(airmiles),end=c(2001,8)))
plot(window(diff(log(airmiles),1),end=c(2001,8)))
acf(window(diff(log(airmiles),1),end=c(2001,8)),lag.max=48)
acf(window(diff(diff(log(airmiles),1),lag=12),end=c(2001,8)),lag.max=48)
pacf(window(diff(diff(log(airmiles),1),lag=12),end=c(2001,8)),lag.max=48)
eacf(window(diff(diff(log(airmiles),1),lag=12),end=c(2001,8)))


arima1=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,0)))
arima1
tsdiag(arima1)
arima2=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)))
arima2
tsdiag(arima2)

arima3=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)),
xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)))
arima3

# half-life
-log(2)/log(.89)

tsdiag(arima3)
detectAO(arima3)
detectIO(arima3)
arima3a=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)),
xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),io=list(25))
arima3a
tsdiag(arima3a)
detectAO(arima3a)
detectIO(arima3a)
arima3b=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)),
xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),io=list(25,84))
arima3b
tsdiag(arima3b)
detectAO(arima3b)
detectIO(arima3b)
arima3c=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)),
xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),io=list(25,84))

arima3d=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1)),
xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),io=list(84), 
xreg=data.frame(AO.25=1*(seq(airmiles)==25)))
arima3d
tsdiag(arima3d)

detectAO(arima3d)
detectIO(arima3d)

air.m3e=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1),
period=12),xtransf=data.frame(I911=1*(seq(airmiles)==69),
I911=1*(seq(airmiles)==69)),
transfer=list(c(0,0),c(1,0)),xreg=data.frame(Dec96=1*(seq(airmiles)==12),
Jan97=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84)),method='ML',
io=list(25,81))
detectAO(air.m3e)
detectIO(air.m3e)
tsdiag(air.m3e)

air.m3e
tsdiag(air.m3e)