simulate ARIMA(1,0,1) and ARIMA(1,1,1) in fable

[Amarjit_Chandhial]

Hi,

I have simulated the following ARIMA(1,0,1) model in fable

y(t) = 3 + 0.7*y(t-1) + e(t) + 0.6*e(t-1)

which I think is correct? Is there a better way?

How do I simulate an ARIMA(1,1,1) model with the SAME parameters, and how do I write down the equation?

library(fpp3)
library(patchwork)

set.seed(1)


nburn <- 100
nobs <- 200
T <- nburn+nobs 

# Define the parameters for the ARIMA(1,0,1) model
phi <- +0.7
theta <- +0.6
constant <- 3

mean <- constant / (1 - phi)
mean

y <- vector()
y[1] <- 0
eps <- rnorm(nburn+nobs)
for(i in 2:T)
  y[i] <- constant + phi * (y[i-1]) + eps[i] + theta * eps[i-1]
sim <- tsibble(idx = seq_len(T), y = y, index = idx)
sim


p1 <- sim |>
  slice(101:n()) |>
  autoplot(y) +
  ggplot2::ggtitle("Simulated Time Plot\ny(t) = c + phi*y(t-1) + e(t) + theta*e(t-1)")


p2 <- sim |>
  slice(101:n()) |>
  ACF(y) |>
  autoplot() +
  ggplot2::ggtitle("ACF\ny(t) = c + phi*y(t-1) + e(t) + theta*e(t-1)")

p1 + p2


# fit ARIMA
fit <- sim |> slice((nburn+1):n()) |>
  model('sim'= ARIMA( y ~ 1 + pdq(1, 0, 1) + PDQ(0, 0, 0), 
                      method = "CSS-ML", optim.control = list(maxit = 1000))
  )	
fit
report(fit)

fit$sim[[1]][[1]][[1]] # constant
fit$sim[[1]][[1]][[6]] # mean

Amarjit

[Amarjit_Chandhial]

Ahhh if I integrate-up, and fit again, taking into account the difference, I get the same coeffs for ar1 and ma1 in both ARIMA(1,0,1) and ARIMA(1,1,1).

# fit ARIMA
fit <- sim |> slice((nburn+2):n()) |>
  model('sim'= ARIMA( y ~ 1 + pdq(1, 0, 1) + PDQ(0, 0, 0), 
                      method = "CSS-ML", optim.control = list(maxit = 1000))
  )	
fit
report(fit)

fit$sim[[1]][[1]][[1]] # constant
fit$sim[[1]][[1]][[6]] # mean

> report(fit)
Series: y 
Model: ARIMA(1,0,1) w/ mean 

Coefficients:
         ar1     ma1  constant
      0.6758  0.5484    3.2227
s.e.  0.0604  0.0788    0.1075

sigma^2 estimated as 0.9973:  log likelihood=-281.38
AIC=570.77   AICc=570.97   BIC=583.94
> 
> fit$sim[[1]][[1]][[1]] # constant
# A tibble: 3 × 5
  term     estimate std.error statistic  p.value
  <chr>       <dbl>     <dbl>     <dbl>    <dbl>
1 ar1         0.676    0.0604     11.2  7.60e-23
2 ma1         0.548    0.0788      6.96 4.73e-11
3 constant    3.22     0.107      30.0  9.69e-76
> fit$sim[[1]][[1]][[6]] # mean

Call:
wrap_arima(x = y, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = period), 
    xreg = xreg, include.mean = FALSE, fixed = fixed, method = method, optim.control = list(
        maxit = 1000))

Coefficients:
         ar1     ma1  constant
      0.6758  0.5484    9.9397
s.e.  0.0604  0.0788    0.3315

sigma^2 estimated as 0.9973:  log likelihood = -281.38,  aic = 570.77

sim <- sim |> mutate(group_by(sim,y), csum=cumsum(y))

# fit ARIMA
fit <- sim |> slice((nburn+1):n()) |>
  model('sim'= ARIMA( csum ~ 1 + pdq(1, 1, 1) + PDQ(0, 0, 0), 
                      method = "CSS-ML", optim.control = list(maxit = 1000))
  )	
fit
report(fit)

# A mable: 1 x 1
                      sim
                  <model>
1 <ARIMA(1,1,1) w/ drift>
> report(fit)
Series: csum 
Model: ARIMA(1,1,1) w/ drift 

Coefficients:
         ar1     ma1  constant
      0.6758  0.5484    3.2228
s.e.  0.0604  0.0788    0.1075

sigma^2 estimated as 1.003:  log likelihood=-281.39
AIC=570.77   AICc=570.98   BIC=583.94

[Amarjit_Chandhial]

Here's 2 more examples.

Simulate ARIMA(1,0,0) with constant

#===============================
library(gratis)
library(fpp3)


#===============================
# Simulate ARIMA(1,0,0) with constant
# -----------------------------------

# Y(t) = 150 + 0.7*Y(t-1) + eps(t)
# --------------------------------

constant <- 150.0
phi1 <- 0.7
mean <- (constant / (1.0 - phi1))
mean # unconditional mean (or isolating behaviour) as you look at the graph


#===============================
# Simulation
ar1 <- gratis::arima_model(p = 1, d = 0, q = 0,
  phi = phi1,
  constant = mean,
  sigma = 1
)

ar1_d0 <- gratis::generate(ar1, length = 200, nseries = 1)
ar1_d0


ar1_d0 |>
  autoplot(value) +
  labs(title="ARIMA(1,0,0)+c\nY(t) = 150 + 0.7*Y(t-1) + eps(t)")

tb_m <- ar1_d0
mean(tb_m$value)


# fit ARIMA (as the model is written down)
fit <- ar1_d0 |>
  model('ar1_d0_sim'= ARIMA( value ~ 1 + pdq(1, 0, 0),
                      fixed = list(ar1 = NA, constant = NA), 
                      transform.pars = FALSE,
                      method = "CSS-ML", optim.control = list(maxit = 1000))
  )	
fit
fit$ar1_d0_sim[[1]] # BOTH constant and mean
report(fit)

# estimated unconditional mean
mean_est <- (((fit[[1]][[1]]$fit$par$estimate[2]) / (1.0 - fit[[1]][[1]]$fit$par$estimate[1]))) 
mean_est


# Base R
sim.ts <- as.ts(ar1_d0)
sim.ts |>
  arima( order=c(1,0,0), include.mean = TRUE,
         fixed = list(ar1 = NA, intercept = NA), 
         transform.pars = FALSE,
         method = "CSS-ML", optim.control = list(maxit = 1000)
  )


# library(forecast)
sim.ts |>
  forecast::Arima( order=c(1,0,0), include.constant = TRUE,
                   fixed = list(ar1 = NA, mean = NA), 
                   transform.pars = FALSE,
                   method = "CSS-ML", optim.control = list(maxit = 1000)
  )

Simulate ARIMA(1,1,0) with constant

#===============================
library(gratis)
library(fpp3)


#===============================
# Simulate ARIMA(1,1,0)
# ---------------------

# d(Y(t)) = 150 + 0.7*d(Y(t-1)) + eps(t) See Brockwell Davis (2016) Introduction to Time Series and Forecasting 3rdEdn, p.158
# --------------------------------------

constant <- 150.0
phi1 <- 0.7
mean <- (constant / (1.0 - phi1))
mean # uncondtional mean (or isolating behaviour) as you look at the graph


#===============================
# Simulation
ari1 <- gratis::arima_model(p = 1, d = 1, q = 0,
  phi = phi1,
  constant = 0, # constant = 0
  sigma = 1
)

ari1_d1 <- gratis::generate(ari1, length = 200, nseries = 1) 
ari1_d1


ari1_d1 |>
  dplyr::mutate(value = value + mean) |> # mean is added onto value
  autoplot(value) +
  labs(title="ARIMA(1,1,0)+c\nd(Y(t)) = 150 + 0.7*d(Y(t-1)) + eps(t)\n(1−0.7*B)*y(t) = 150 + e(t), y(t) = Y(t)-Y(t-1)")

tb_m <- ari1_d1 + mean # mean is added onto value
mean(tb_m$value)


# fit ARIMA (as the model is written down)
fit <- ari1_d1 |>
  model('ari1_d1_sim'= ARIMA( value ~ 0 + pdq(1, 1, 0),
                              fixed = list(ar1 = NA), 
                              transform.pars = FALSE,
                              method = "CSS-ML", optim.control = list(maxit = 1000))
  )	
fit
fit$ari1_d1_sim[[1]] # BOTH constant and mean
report(fit)


# Base R
sim.ts <- as.ts(ari1_d1)
b.fit <- sim.ts |>
  stats::arima( order=c(1,1,0),
                fixed = list(ar1 = NA), 
                transform.pars = FALSE,
                method = "CSS-ML", optim.control = list(maxit = 1000)
  )
b.fit


# library(forecast)
f.fit <- sim.ts |>
  forecast::Arima( order=c(1,1,0),
                   fixed = list(ar1 = NA), 
                   transform.pars = FALSE,
                   method = "CSS-ML", optim.control = list(maxit = 1000)
  )
f.fit

The questions being

• Are they simulated correctly?
• Do I add on the mean after the simulation, therefore the constant is not estimated for ARIMA(1,1,0) with constant?
• Is the ARIMA(1,1,0) with constant written down correctly?