1 Simulate data
First, we simulate the innovation distribution. Note that, for demonstration purposes, we choose a small sample size. Ideally, the sample size should be larger to capture GARCH effects.
## Simulate innovations
n <- 200 # sample size
d <- 2 # dimension
nu <- 3 # degrees of freedom for t
tau <- 0.5 # Kendall's tau
th <- iTau(ellipCopula("t", df = nu), tau) # corresponding parameter
cop <- ellipCopula("t", param = th, dim = d, df = nu) # define copula object
set.seed(271) # reproducibility
U <- rCopula(n, cop) # sample the copula
nu. <- 3.5 # degrees of freedom for the t margins
Z <- sqrt((nu.-2)/nu.) * qt(U, df = nu.) # margins must have mean 0 and variance 1 for ugarchpath()!Now we simulate two ARMA(1,1)-GARCH(1,1) processes with these copula-dependent innovations. To this end, recall that an ARMA(\(p_1\),\(q_1\))-GARCH(\(p_2\),\(q_2\)) model is given by \[\begin{align} X_t &= \mu_t + \epsilon_t\ \text{for}\ \epsilon_t = \sigma_t Z_t,\\ \mu_t &= \mu + \sum_{k=1}^{p_1} \phi_k (X_{t-k}-\mu) + \sum_{k=1}^{q_1} \theta_k \epsilon_{t-k},\\ \sigma_t^2 &= \alpha_0 + \sum_{k=1}^{p_2} \alpha_k (X_{t-k}-\mu_{t-k})^2 + \sum_{k=1}^{q_2} \beta_k \sigma_{t-k}^2. \end{align}\]
## Fix parameters for the marginal models
fixed.p <- list(mu = 1,
ar1 = 0.5,
ma1 = 0.3,
omega = 2, # alpha_0 (conditional variance intercept)
alpha1 = 0.4,
beta1 = 0.2)
meanModel <- list(armaOrder = c(1,1))
varModel <- list(model = "sGARCH", garchOrder = c(1,1)) # standard GARCH
uspec <- ugarchspec(varModel, mean.model = meanModel,
fixed.pars = fixed.p) # conditional innovation density (or use, e.g., "std")
## Simulate ARMA-GARCH models using the dependent innovations
## Note: ugarchpath(): simulate from a spec; ugarchsim(): simulate from a fitted object
X <- ugarchpath(uspec,
n.sim = n, # simulated path length
m.sim = d, # number of paths to simulate
custom.dist = list(name = "sample", distfit = Z)) # passing (n, d)-matrix of *standardized* innovations
## Extract the resulting series
X. <- fitted(X) # X_t = mu_t + eps_t (simulated process)
sig.X <- sigma(X) # sigma_t (conditional standard deviations)
eps.X <- X@path$residSim # epsilon_t = sigma_t * Z_t (residuals)
## Basic sanity checks :
stopifnot(all.equal(X., X@path$seriesSim, check.attributes = FALSE),
all.equal(sig.X, X@path$sigmaSim, check.attributes = FALSE),
all.equal(eps.X, sig.X * Z, check.attributes = FALSE))
## Plot (X_t) for each margin
matplot(X., type = "l", xlab = "t", ylab = expression(X[t]~"for each margin"))
