phaser - 5 months ago 26

R Question

I am working with

`arima0()`

`co2`

`arima0()`

`fitted()`

`curve()`

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)
```

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)
```

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)
```

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)
```

Oh, this is just gorgeous!