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

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.