I am studying $ARFIMA$ processes, and I have some trouble understanding how I can simulate it in R. I know that there are functions and packages to do so, but I'd like to understand the mathematical way of doing so.
If I have for instance an $ARMA(p,q)$ and a sample of size $T$, it follows that for each $t \in \{1, 2, \ldots T\}$:
$$ X_t = \sum_{j=1}^{p} \phi_j X_{t-j} + \epsilon_t + \sum_{i=1}^{q} \theta_i \epsilon_{t-i} $$ $\epsilon \sim WN(0, \sigma^2)$, obeying all the stationarity and invertibility conditions, meaning that the roots of the polinomials $\Phi(B) = (1-\phi_1 B - \ldots \phi_p B^p) $ and $ \Theta(B) = (1 + B\theta_1 + B^2 \theta_2 + \ldots + B^q \theta_q)$ are not in the unit circle and $B$ is the shift operator, meaning: $B^p X_t = X_{t-p}$.
So if I want to simulate for instance an $ARMA(1,1)$, I'd do the following:
- Simulate an gaussian white noise $E_t \sim N(,\sigma_E^2) $
- Set a $X_0=0$ for instance (this won't matter lately).
- Evaluate $X_t = \phi X_{t-1} + \theta E_t $, for given parameters $(\phi, \theta, \sigma^2_E)$.
- If I want my sample to be of size $T$, I'd simulate $2T$ values and discard the first $T$, to avoind being biased.
Now $ARFIMA$ models are such that:
$$ \Phi(B)(1-B)^d X_t = \Theta(B) e_t $$, where $\Phi(B)$ and $\Theta(B)$ are the above defined polynomials and $B$ is the shift operator. The interesting thing is $d$ as a number such that $d \in (-0.5,0.5)$. How could I simulate an ARFIMA process as in the routine I above defined? Because my routine is so simple I coudl reproduce it even in an Excell sheet.
Thanks!