Solved – Simulation and mathematical notation for ARIMA(0,1,1) with drift

arimaforecastingrsimulationtime series

I am attempting to write the mathematical model for and also simulate an MA(1) process that has drift (in R).

I have referenced ARIMA (0,1,1) or (0,1,0) – or something else?, Simulation of forecasted values in ARIMA (0,1,1), and Fitting ARIMA with a drift on R. I am aware that the "forecast" package that I am using addresses the issues relating to mean reporting here: http://www.stat.pitt.edu/stoffer/tsa2/Rissues.htm.

My understanding (e.g. from example bottom of page on the very useful https://www.otexts.org/fpp/8/7 page on the otext site) is that an MA(1) model (on a series that for which the first difference is stationary) can be written as:
$$
(1-B)^d(y_t- \mu t) = (1 – \theta B) e_t
$$
where $B$ is the backshift operator, $d$ is the order of differencing (here is 1), $y_t$ is the original series indexed by time $t$, $\theta$ is the MA1 parameter and $e_t$ is the error indexed at time $t$.

From the above otext site, I note:

the inclusion of a constant in a non-stationary ARIMA model is equivalent to inducing a polynomial trend of order $d$ in the forecast function.

So while my forecasts are correct, Q.1. How do I represent/differentiate this single order trend in the equation from the equation for an MA(1) model without drift? I ask this because expanding the above formula just produces the exact same equation/model I expect for an MA(1) without drift, i.e. (and apologies in advance if there is just something obviously wrong with my maths):
$$
\begin{aligned}
y_t – \mu t – y_{t-1} + \mu t &= e_t – \theta e_{t-1} \\
y_t – y_{t-1} &= e_t – \theta e_{t-1} \\
\Delta y_t &= e_t – \theta e_{t-1} \\
\end{aligned}
$$

In terms of simulating an MA(1) with drift I used the following, Q.2. Is this the correct approach or is there a more direct/accurate way to do this?:

library(lattice)
library(forecast)
set.seed(325354)
temp <- arima.sim(n=100, list(order=c(0,1,1), ma=.75), sd = 2.7) 
mean(temp)

ts.sim.orig <- ts(temp, start=c(2011,1), frequency=12)
xyplot(ts.sim.orig)
mean(ts.sim.orig)
# Stationary
kpss.test(ts.sim.orig)

# Drift term being 0.3
ts.sim <- ts(1:length(ts.sim.orig) * 0.3 + temp, start=c(2011,1), frequency=12)
xyplot(ts.sim)
mean(ts.sim)
# Not stationary
kpss.test(ts.sim)
# Stationary
kpss.test(diff(ts.sim, 1))
# Non zero mean on first differnce
mean(diff(ts.sim, 1))
# Plot is reasonable?
xyplot(diff(ts.sim, 1))

# First fit without drift:
summary(fit99 <- Arima(ts.sim, order=c(0,1,1), include.constant = F)) 
# Forecast is flatline as expected. Great.
plot(forecast(fit99))

# The drift is close to the mean from: mean(diff(ts.sim, 1))
summary(fit99 <- Arima(ts.sim, order=c(0,1,1), include.constant = T)) 
# Forecast includes trend
plot(forecast(fit99))

y <- ts.sim
x <- 1:length(ts.sim)
# Roughly equal to drift term as expected
summary(lm(y ~ x))

# Aside, this don't seem to recover the drift parameter too well...
summary(fit99 <- auto.arima(ts.sim, d = 1, trace = T, stepwise=FALSE, approximation=FALSE, seasonal = F))
plot(forecast(fit99))

Finally, does the approach to the simulation above answer my first question? That is, is the appropriate representation of the MA(1) with drift (and presumably assuming a zero intercept) equal to the following equation?
$$
\Delta y_t = \beta t + e_t – \theta e_{t-1}
$$
Thank you.

Best Answer

Regarding Q1 : The difference operator applied to µt gives µt – µ(t-1) = µ. Hence the drift parameter µ will appear in your equation. Maybe better to choose another symbol for the drift parameter, since µ is traditionally used for the mean of the time series.

Related Question