library(fable)
library(tsibble)
library(dplyr)

#create the variable to store ozone values 
#Box & Jenkins Series R
y = c(2.63, 1.94, 3.38, 4.92, 6.29, 5.58, 5.5, 4.71, 6.04, 
7.13, 7.79, 3.83, 3.83, 4.25, 5.29, 3.75, 4.67, 5.42, 6.04, 5.71, 
8.13, 4.88, 5.42, 5.5, 3, 3.42, 4.5, 4.25, 4, 5.33, 5.79, 6.58, 
7.29, 5.04, 5.04, 4.48, 3.33, 2.88, 2.5, 3.83, 4.17, 4.42, 4.25, 
4.08, 4.88, 4.54, 4.25, 4.21, 2.75, 2.42, 4.5, 5.21, 4, 7.54, 
7.38, 5.96, 5.08, 5.46, 4.79, 2.67, 1.71, 1.92, 3.38, 3.98, 4.63, 
4.88, 5.17, 4.83, 5.29, 3.71, 2.46, 2.17, 2.15, 2.44, 2.54, 3.25, 
2.81, 4.21, 4.13, 4.17, 3.75, 3.83, 2.42, 2.17, 2.33, 2, 2.13, 
4.46, 3.17, 3.25, 4.08, 5.42, 4.5, 4.88, 2.83, 2.75, 1.63, 3.04, 
2.58, 2.92, 3.29, 3.71, 4.88, 4.63, 4.83, 3.42, 2.38, 2.33, 1.5, 
2.25, 2.63, 2.96, 3.46, 4.33, 5.42, 4.79, 4.38, 4.54, 2.04, 1.33, 
2.04, 2.81, 2.67, 4.08, 3.9, 3.96, 4.5, 5.58, 4.52, 5.88, 3.67, 
1.79, 1.71, 1.92, 3.58, 4.4, 3.79, 5.52, 5.5, 5, 5.48, 4.81, 
2.42, 1.46, 1.71, 2.46, 2.42, 1.79, 3.63, 3.54, 4.88, 4.96, 3.63, 
5.46, 3.08, 1.75, 2.13, 2.58, 2.75, 3.15, 3.46, 3.33, 4.67, 4.13, 
4.73, 3.42, 3.08, 1.79, 1.96, 1.63, 2.75, 3.06, 4.31, 3.31, 3.71, 
5.25, 3.67, 3.1, 2.25, 2.29, 1.25, 2.25, 2.67, 3.23, 3.58, 3.04, 
3.75, 4.54, 4.46, 2.83, 1.63, 1.17, 1.79, 1.92, 2.25, 2.96, 2.38, 
3.38, 3.38, 3.21, 2.58, 2.42, 1.58, 1.21, 1.42, 1.96, 3.04, 2.92, 
3.58, 3.33, 4.04, 3.92, 3.08, 2, 1.58, 1.21)

#Create the intervention dummy
x = 1:216

for (i in 1:228){
    if (i<=60){
        x[i] = 0
    }
    else x[i] = 1
}

#create the tsibble object for model fitting
dt = tsibble(date=rep(yearmonth("1955 Jan") + 0:215), index=date)
extra = tsibble(date=rep(yearmonth("1973 Jan") + 0:11), index=date)

dt$ozone = y
dt$x1 = x[1:216]

#create the tsibble object for future x values
extra$x1 = x[217:228]

#define the ARIMAX model
fit = dt %>% model(arima=ARIMA(ozone~ 0 + x1 + pdq(0,0,1) + PDQ(0,1,1), method="ML"))

#report parameter estimates
fit %>% report()

#report forecasts
fit %>% forecast(new_data=extra)