
/*************** Simulating ARCH variates  *************************
Author: Mico Loretan  loretanm@frb.gov
Last Update: 14 Oct 96            Previous: 19 Aug 1996
This code is provided without guarantee for public non-commercial use.

Contains three procs for time series simulation:
  "arch()", "garch()" and "igarch()".
***********************************************************************/

proc(2)= arch(n,s2,phi);
    /*
    :: Mico Loretan (loretanm@frb.gov)
    :: August, 1996
    ::
    :: Proc to generate ARCH(q) type innovations with finite
    :: unconditional variance.  Proc does not check explicitly
    :: for sufficient conditions to guarantee positive cond.
    :: variances.
    ::
    :: The ARCH(q) process is:
    ::    e_t ~ N(0,h_t^2), where
    ::    h_t^2 = s2*(1-sumc(phi)) +
    ::              \sum_{j=1}^q \phi_j e_{t-j}^2.
    ::
    :: Inputs:    n     number of observations to be generated
    ::            s2    unconditional variance
    ::            phi   coefficient vector (q x 1)
    ::
    :: Outputs:   e     n x 1 vector of ARCH(q) innovations
    ::            h2    associated conditional variance process
    ::                  (n x 1)
    ::
    :: Remark:    We generate ini=2q "additional" observations,
    ::            which are discarded as startup/transitional
    ::            values.  This is ad-hoc; if deemed necessary,
    ::            increase "ini".
    */

    local q,ini, a0, h2,e, t;

    /* first, check for valid inputs */

    if ( (sumc(phi)>=1) or (s2<=0) );
    errorlog "Unable to generate ARCH(q) process with positive variance";
    end;
    endif;

    q   = rows(phi);    /* order of ARCH process */
    ini = 2*q;          /* increase this if desired; do not
                           change the value of ini below */

    phi = rev(phi);         /* reverse order of coefficients */
    a0  = (1-sumc(phi))*s2;

    h2 = zeros(n+ini,1);    /* initialize output vectors */
    e  = zeros(n+ini,1);

    e[1:pq]  = sqrt(s2)*rndn(q,1); /* need q initial obs. */
    h2[1:pq] = ones(q,1)*s2;       /* for first q values, h2:=s2*/

    t = q+1;
    do until (t>n+ini);
        h2[t] = a0 + phi'((e[t-q:t-1])^2);
        e[t]  = sqrt(h2[t]) * rndn(1,1);
        t     = t+1;
    endo;

    e  = e[ini+1:ini+n];    /* discard ini startup values */
    h2 = h2[ini+1:ini+n];

    retp(e,h2);
endp; /* arch() */

proc(2)= garch(n,s2,psi,phi);
    /*
    :: Mico Loretan (loretanm@frb.gov)
    :: August 1996
    ::
    :: Proc to generate GARCH(p,q) type innovations with finite
    :: unconditional variance.  Proc does not check explicitly
    :: for sufficient conditions to guarantee positive cond'l
    :: variances.  Furthermore, the proc does not check that at
    :: least one of the elements of "phi" is nonzero.
    ::
    :: The GARCH(p,q) model is:
    ::      e_t ~ N(0,h_t^2), where
    ::      h_t^2 = s2*(1-sumc(psi)-sumc(phi)) +
    ::              \sum_{i=1}^p \psi_i h_{t-i} +
    ::              \sum_{j=1}^q \phi_j e_{t-j}^2.
    ::
    :: Inputs:    n     number of observations to be generated
    ::            s2    unconditional variance
    ::            psi   coefficient vector for lagged conditional
    ::                  variances (p x 1)
    ::            phi   coefficient vector for lagged squared
    ::                  innovations (q x 1)
    :: Outputs:   e     vector of GARCH(p,q) innovations (n x 1)
    ::            h2    associated conditional variance process
    ::                  (n x 1)
    :: Remark:    We generate ini=2*maxc(p|q) "additional"
    ::            observations, which are discarded as startup/
    ::            transitional values.  This is ad-hoc; if deemed
    ::            necessary, increase the value of the variable
    ::            "ini".
    */
    local p,q,pq,ini, a0, h2,e, t;

    /* first, check for valid inputs */

    if ( (sumc(phi)+sumc(psi) >= 1) or (s2 <= 0) );
    errorlog "Unable to generate GARCH(p,q) process with positive variance";
    end;
    endif;

    p   = rows(psi);     /* order of GARCH(p,q) process */
    q   = rows(phi);
    pq  = maxc(p|q);
    ini = 2*pq;          /* increase this if desired; do not
                            change the value of ini below */

    psi = rev(psi);      /* reverse order of coefficients */
    phi = rev(phi);
    a0  = (1-sumc(psi)-sumc(phi)) * s2;

    h2 = zeros(n+ini,1); /* initialize output vectors */
    e  = zeros(n+ini,1);

    e[1:pq]  = sqrt(s2)*rndn(pq,1); /* need pq initial obs. */
    h2[1:pq] = ones(pq,1)*s2;       /* for first pq values,h2:=s2 */

    t = pq+1;
    do until (t>n+ini);
        h2[t] = a0 + psi'(h2[t-p:t-1]) + phi'(e[t-q:t-1])^2;
        e[t]  = sqrt(h2[t]) * rndn(1,1);
        t     = t+1;
    endo;

    e  = e[ini+1:ini+n];    /* discard ini startup values */
    h2 = h2[ini+1:ini+n];

    retp(e,h2);
endp; /* garch() */





proc(2)= igarch(n,a0,psi,phi);
    /*
    :: Mico Loretan (loretanm@frb.gov)
    :: August, 1996
    ::
    :: Proc to generate IGARCH(p,q) type innovations.
    ::
    :: The IGARCH(p,q) model is:
    ::      e_t ~ N(0,h_t^2), where
    ::      h_t^2 = a0 + \sum_{i=1}^p \psi_i h_{t-i} +
    ::              \sum_{j=1}^q \phi_j e_{t-j}^2
    ::      and \sum \psi_i + \sum \phi_j = 1.
    ::
    :: Inputs:    n     number of observations to be generated
    ::            a0    constant, > 0
    ::            psi   coefficient vector for lagged conditional
    ::                  variances (p x 1)
    ::            phi   coefficient vector for lagged squared
    ::                  innovations (q x 1)
    :: Outputs:   e     vector of GARCH(p,q) innovations (n x 1)
    ::            h2    associated conditional variance process
    ::                  (n x 1)
    :: Remark:    We generate ini=2*maxc(p|q) "additional"
    ::            observations, which are discarded as startup/
    ::            transitional values.  This is ad-hoc; if deemed
    ::            necessary, increase the value of the variable
    ::            "ini".
    */
    local p,q, pq,ini, h2,e, t;

    /* first, check for valid inputs */

    if ( (sumc(phi)+sumc(psi) ne 1) or (a0 le 0) );
    errorlog "Unable to generate proper IGARCH(p,q) process";
    end;
    endif;

    p   = rows(psi);     /* order of IGARCH(p,q) process */
    q   = rows(phi);
    pq  = maxc(p|q);
    ini = 4*pq;          /* increase this if desired; do not
                            change the value of ini below */

    psi = rev(psi);      /* reverse order of coefficients */
    phi = rev(phi);

    h2 = zeros(n+ini,1); /* initialize output vectors */
    e  = zeros(n+ini,1);

    e[1:pq]  = sqrt(a0)*rndn(pq,1); /* need pq initial obs. */
    h2[1:pq] = ones(pq,1)*a0; /* for first pq values, h2:=a0 */

    t = pq+1;
    do until (t>n+ini);
        h2[t] = a0 + psi'(h2[t-p:t-1]) + phi'(e[t-q:t-1])^2;
        e[t]  = sqrt(h2[t]) * rndn(1,1);
        t     = t+1;
    endo;

    e  = e[ini+1:ini+n];    /* discard ini startup values */
    h2 = h2[ini+1:ini+n];

    retp(e,h2);
endp; /* igarch() */
