phaser phaser - 1 month ago 12
R Question

ARIMA modelling, prediction and plotting with CO2 dataset in R

I am working with

arima0()
and
co2
. I would like to plot
arima0()
model over my data. I have tried
fitted()
and
curve()
with no success.

Here is my code:

###### Time Series

# format: time series
data(co2)

# format: matrix
dmn <- list(month.abb, unique(floor(time(co2))))
co2.m <- matrix(co2, 12, dimnames = dmn)

co2.dt <- pracma::detrend(co2.m, tt = 'linear')
co2.dt <- ts(as.numeric(co2.dt), start = c(1959,1), frequency=12)

# first diff
co2.dt.dif <- diff(co2.dt,lag = 12)

# Second diff
co2.dt.dif2 <- diff(co2.dt.dif,lag = 1)


With the data prepared, I ran the following
arima0
:

results <- arima0(co2.dt.dif2, order = c(2,0,0), method = "ML")
resultspredict <- predict(results, n.ahead = 36)


I would like to plot the model and the prediction. I am hoping there is a way to do this in base R. I would also like to be able to plot the predictions as well.

Answer

Session 2: A better idea for time series modelling

In previous session, we have seen the difficulty in prediction if we separate de-trending and modelling of de-trended series. Now, we try to combine those two stages in one go.

The seasonal pattern of co2 is really strong, so we need a seasonal differencing anyway:

data(co2)
co2dt <- diff(co2, lag = 12)
par(mfrow = c(1,2)); ts.plot(co2dt); acf(co2dt)

enter image description here

After this seasonal differencing, co2dt does not look stationary. So we need further a non-seasonal differencing.

co2dt.dif <- diff(co2dt)
par(mfrow = c(1,2)); ts.plot(co2dt.dif); acf(co2dt.dif)

enter image description here

The negative spikes within season and between season suggest that a MA process is needed for both. I will not work with co2dt.dif; we can work with co2 directly:

fit <- arima0(co2, order = c(0,1,1), seasonal = c(0,1,1))
acf(fit$residuals)

enter image description here

Now the residuals are perfectly uncorrelated! So we have an ARIMA(0,1,1)(0,1,1)[12] model for co2 series.


As usual, fitted values are obtained by subtracting residuals from data:

co2fitted <- co2 - fit$residuals

Predictions are made by a single call to predict:

co2pred <- predict(fit, n.ahead = 36, se.fit = FALSE)

Let's plot them together:

ts.plot(co2, co2fitted, co2pred, col = 1:3)

enter image description here

Oh, this is just gorgeous!