Seasonal ARMA and ARIMA Models

 Data: Value of fresh fruit imports into Canada (millions $), collected
monthly.

> par(mfrow=c(3,1))
> tsplot(fruitimp.ts)
> acf(fruitimp.ts)
> acf(fruitimp.ts, type = "par")

 What do these plots suggest?

 # First difference to remove linear trend

> fruitimp.ts.diff <- diff(fruitimp.ts, lag = 1, differences = 1)    

> par(mfrow=c(2,2))
> acf(fruitimp.ts.diff, lag.max = 48)
> acf(fruitimp.ts.diff, type = "par", lag.max = 48)

# Although ACF in seasonal component is going down, it's not decaying
# quickly

# Nonstationary in seasonal component?
# Try first difference in months (seasons). 

> fruitimp.ts.diff.seas <- diff(fruitimp.ts.diff,lag = 12,differences=1) 
> acf(fruitimp.ts.diff.seas, lag.max = 48)
> acf(fruitimp.ts.diff.seas, type = "par", lag.max = 48)

 # Suitable ARMA model after adjusting for nonstationarity?



################################################

 # MA(1)

> arima.mle(fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas),
   model = list(ma = .2))

Call: arima.mle(x = fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas),
model = list(  ma = 0.2)) 
Method:  Maximum Likelihood 
Model :  0 0 1 

Coefficients:
        MA : 0.87365

Variance-Covariance Matrix:
             ma(1) 
ma(1) 0.0009431436

Optimizer has  converged
Convergence Type: relative function convergence
AIC: 5243.7511

# AR(2)

> arima.mle(fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas),
 model = list(ar = c(.2, .2)))

Call: arima.mle(x = fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas),
model = list(  ar = c(0.2, 0.2))) 
Method:  Maximum Likelihood 
Model :  2 0 0 

Coefficients:
        AR : -0.69216 -0.38236

Variance-Covariance Matrix:
            ar(1)       ar(2) 
ar(1) 0.003428918 0.001716879
ar(2) 0.001716879 0.003428918

Optimizer has  converged
Convergence Type: relative function convergence
AIC: 5229.20259

######################################################

 Try (0,1) x (0,1)_12, to allow for MA(1) in seasonal component:
 Not specifying starting values: Splus will choose 0 as default

> arima.mle(fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas),
  model = list(list(order = c(0,0,1)),
  list(order = c(0,0,1), period=12)))

Call: arima.mle(x = fruitimp.ts.diff.seas - mean(fruitimp.ts.diff.seas), 
model = list(list(order = c(0, 0, 1)), list(order = c(0, 0, 1), period = 12)))
Method:  Maximum Likelihood
Model :
 
Coefficients: WHERE ARE THEY? SEE EXAMPLE BELOW TO SEE HOW TO FIND THEM.
 
Variance-Covariance Matrix:
               ma(1)        ma(12)
 ma(1)  0.0004867111 -0.0003614092
ma(12) -0.0003614092  0.0032034370
 
Optimizer has  converged
Convergence Type: relative function convergence
AIC: 5197.24426
#########################################

 # Fitting (0,1,1) x (0,1,1)_12 directly to data. This is equivalent to
 # what is done above (fitting (0,1) x (0,1)_12 to differenced data).

> fruit.model <- arima.mle(fruitimp.ts - mean(fruitimp.ts),
 model = list(list(order = c(0,1,1)), list(order = c(0,1,1),period=12)))

> fruit.model # Some output deleted
Model :
 
Coefficients:
 
Variance-Covariance Matrix:
               ma(1)        ma(12)
 ma(1)  0.0003554065 -0.0003239262
ma(12) -0.0003239262  0.0032254111
 
Convergence Type: relative function convergence
AIC: 5196.49708

> fruit.model$model # To see the estimated coefficients
[[1]]:
[[1]]$order:
[1] 0 1 1
 
[[1]]$ndiff:
[1] 1
 
[[1]]$ma:
[1] 0.9586232    # Close to non-invertible

[[2]]:
[[2]]$order:
[1] 0 1 1
 
[[2]]$period:
[1] 12
 
[[2]]$ndiff:
[1] 1
 
[[2]]$ma:
[1] 0.5143206

###########################################################

 Forecasting for next 12 months in future:

> fruit.fore <- arima.forecast(fruitimp.ts - mean(fruitimp.ts),
  n=12, model = fruit.model$model)

> fruit.fore
$mean:          # Forecasts
 [1]  30927.98  34119.87  57128.69  54528.04  80809.81 101787.61  84297.46
 [8]  70456.03  52166.08  36006.46  58948.29  51334.52
 
$std.err:      # Standard Errors of Forecasts
 [1] 7411.991 7418.333 7424.669 7431.001 7437.327 7443.647 7449.962 7456.272
 [9] 7462.577 7468.876 7475.170 7481.458

 # To get forecasts on original scale: add mean back in:

> fruit.fore$mean + mean(fruitimp.ts)
 [1] 126957.2 130149.1 153157.9 150557.2 176839.0 197816.8 180326.6 166485.2
 [9] 148195.3 132035.6 154977.5 147363.7