"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/addons/gig/gig_mle.inp" (11 Nov 2020, 13454 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Charmm source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "gig_mle.inp": 2020d_vs_2020e.

    1 function series ln_pdf_skt(series *u, scalar df, scalar ht_skew)
    2     series ret = NA
    3 
    4     if (df>2) 
    5     # sqrt(pi) = 1.77245385090551602729816748334
    6         scalar q = lngamma((df+1)/2) - lngamma(df/2)
    7     scalar c = exp(q)/(sqrt(df-2)*1.77245385090551602729816748334)
    8 
    9     scalar a = 4 * ht_skew * c * ((df-2)/(df-1))
   10     scalar b = sqrt(1 + 3*ht_skew^2 - a^2)
   11     series d = (b*u + a)
   12     d = (d<0) ? d/(1-ht_skew) : d/(1+ht_skew)
   13     ret = log(b) + log(c) - ((df+1)/2) * log(1+(d^2/(df-2)))
   14     endif
   15 
   16     return ret
   17 end function
   18 
   19 function series ln_pdf_skged(series *x, scalar ni, scalar ta)
   20     scalar p  = 1/ni
   21 
   22     lgp  = lngamma(p)
   23     lg2p = lngamma(2*p)
   24     lg3p = lngamma(3*p)
   25     tap1 = 1 + ta    
   26     tam1 = 1 - ta    
   27 
   28     scalar beta = 0.5 * exp(lg3p - lgp) * (tap1^3 + tam1^3) - \
   29         4*ta^2 * exp(2 * (lg2p - lgp))
   30     beta = sqrt(beta)
   31 
   32     # m: mode
   33     scalar m = - 2*ta/beta * exp( lg2p - lgp )
   34     scalar lnorm = log(0.5 * beta) - lngamma(p+1)
   35 
   36     series z = (x<m) ? (m-x)*beta/tam1 : (x-m)*beta/tap1 
   37     ret = lnorm - (z^ni)
   38     return ret
   39 end function
   40 
   41 
   42 /*
   43  * One-size-fits-all Loglikelihood function
   44  */
   45 
   46 function series gig_loglik(series *e, series *h, scalar distrType, 
   47                matrix addpar, series *de, series *dh, 
   48                matrix *dd)
   49     series ret = NA
   50     de = NA
   51     dh = NA
   52     dd = {}
   53 
   54     if distrType == 0 # Normal
   55     series u = e/h
   56         ret = -.91893853320467274177 - 0.5*(log(h) + e*u)
   57     de = -u
   58     dh = -0.5/h * (1 - e*u)
   59 
   60     elif distrType == 1 # Student's t
   61     scalar ni = addpar[1]
   62     if (ni>2)
   63 
   64             # series hadj = sqrt(h*(1-2/ni))
   65             # series u = e/hadj
   66         # ret = log(pdf(t, ni, u)/hadj)
   67 
   68         series e2 = e*e
   69         scalar K1 = lngamma((ni+1)/2) - lngamma(ni/2) - 0.5*log($pi*(ni-2))
   70         series ret = K1 - 0.5 * log(h) - 0.5*(ni+1) * log(1 + e2/(h*(ni-2)))
   71 
   72         series den = e2 + (ni-2)*h
   73         de = - (ni + 1) * e / den
   74         dh = 0.5/h  * ((ni + 1)* e2 / den - 1)
   75 
   76         scalar k1 = digamma((ni+1)/2) - digamma(ni/2) - 1/(ni-2)
   77         series s1 = (ni + 1)/(ni - 2) * e2/den
   78         series s2 = log(1 + e2/(h*(ni-2)))
   79 
   80         matrix dd = 0.5*(k1 + s1 - s2)
   81     endif
   82 
   83     elif distrType == 2 # GED
   84     scalar ni = addpar[1]
   85     if (ni>0) 
   86         scalar p = 1/ni
   87         scalar lg1 = lngamma(p)
   88         scalar lg3 = lngamma(3*p)
   89         
   90         scalar lC = log(ni/2) + 0.5*(lg3 - 3*lg1)
   91         scalar k = exp(0.5*(lg1-lg3)) * (0.5^p)
   92         series u = abs(e)/(k*sqrt(h))
   93         ret = lC - 0.5*(u^ni + log(h))
   94     endif
   95     
   96     elif distrType == 3 # Skewed-T
   97     scalar ni = addpar[1]
   98     if (ni>2) 
   99         series u = e/sqrt(h)
  100         alpha    = tanh(addpar[2])
  101         ret      = ln_pdf_skt(&u, ni, alpha) - 0.5*log(h)
  102     endif
  103 
  104     elif distrType == 4 # Skewed-GED
  105     scalar ni = addpar[1]
  106     if (ni>0) 
  107         series u = e/sqrt(h)
  108         alpha    = tanh(addpar[2])
  109         ret      = ln_pdf_skged(&u, ni, alpha) - 0.5*log(h)
  110     endif
  111 
  112     endif
  113 
  114     return ret
  115 end function
  116 
  117 /*
  118  * Computes conditional variance and errors for the APARCH model.
  119  * Returns an error code if any of the parameters is out of its domain.
  120  */
  121 
  122 function scalar aparchFilter(series *depVar, series *h, series *e, 
  123                  matrix *mReg, matrix *vX, 
  124                  scalar p, scalar q, matrix *parameters,
  125                  scalar is_asymmetric,
  126                  matrix *deriv_h[null], matrix *deriv_e[null])
  127 
  128    nmX = cols(mReg)
  129    nvX = cols(vX)
  130    scalar base = nmX + nvX
  131 
  132    a0pos = base + 1
  133    a1pos = a0pos + q -1
  134    g0pos = a1pos + 1
  135    g1pos = g0pos + q - 1
  136    b0pos = g1pos + 1
  137    b1pos = b0pos + p - 1
  138    dpos  = b1pos + 1
  139 
  140    matrix avec = {0}
  141    matrix gvec = {0}
  142    matrix bvec = {0}
  143 
  144    matrix omegas = parameters[nmX+1:base]
  145 
  146    if q > 0
  147        matrix avec = parameters[a0pos:a1pos]
  148        matrix gvec = parameters[g0pos:g1pos]
  149    endif
  150 
  151    if p > 0
  152        matrix bvec = parameters[b0pos:b1pos]
  153    endif
  154 
  155    delta = parameters[dpos]
  156 
  157    err = 0
  158 
  159    # Checking
  160    # checks on alpha & beta are disabled
  161    err = err || ((nvX==1) && omegas[1] < 0)
  162 #  err = err || (sumc(avec) + sumc(bvec)) > 1
  163 #  err = err || minc(avec | bvec) < 0
  164    err = err || delta <= 0
  165    # shape? gamma?
  166 
  167    minh = 0
  168 
  169    if err == 0
  170        # handle the conditional mean first
  171        if nmX > 0
  172            series e = depVar - (mReg * parameters[1:nmX])
  173        else
  174            series e = depVar
  175        endif
  176        scalar e0 = mean(e^2)
  177 
  178        matrix tmp_ae = mlag({abs(e)}, seq(1,q), e0 ^ (1/delta))
  179        if is_asymmetric == 1
  180        elag = mlag({e}, seq(1,q))
  181            tmp_ae -= elag .* gvec'
  182        endif
  183 
  184        # Raising a negative value to a non-integer 
  185        # produces a complex number
  186        err = (delta != floor(delta)) && (minc(minr(tmp_ae)) < 0) 
  187        if (!err)
  188        Kd = tmp_ae .^ delta
  189        series h = Kd * avec
  190        
  191        # Var regressors
  192        if (nvX>1)
  193            series h += vX * omegas 
  194        else
  195            series h += omegas[1]
  196        endif
  197 
  198        if p>0
  199            h = filter(h, null, bvec, e0)
  200        endif
  201 
  202        err = min(h)<0 || max(!ok(h))
  203        # the loglikelihood function needs sigma^2
  204        if !err && (delta != 2)
  205            h = h^(2/delta)
  206        endif
  207        endif
  208    endif
  209 
  210    if (!err && exists(deriv_h) && exists(deriv_e))
  211        # FIXME: incomplete and experimental --------------------------------
  212        #
  213        # deriv_e and deriv_h should contain (eventually), the derivatives
  214        # of (doh!) e and h, WITH RESPECT TO THE PARAMETERS
  215        # what we have atm is a rough attempt to have it working in the GARCH 
  216        # case; we'll see about generalising it later
  217 
  218        # ----------- mean eq. --------------
  219        if nmX > 0
  220            deriv_e = -mReg
  221        else
  222            deriv_e = {}
  223        endif
  224 
  225        scalar zcols = nvX + p + q * (1+is_asymmetric)
  226        deriv_e ~= zeros(rows(vX), zcols)
  227 
  228        # ----------- var eq. ---------------
  229        # omega comes 2nd from last 
  230        # mu comes last 
  231        
  232        # alphas
  233        me = Kd
  234        matrix eeff = tmp_ae.^(delta-1)
  235 
  236        if is_asymmetric == 1 # deltas
  237            mh = -delta .* eeff .* ( elag .* avec' )
  238        me ~= mh
  239        endif
  240 
  241        if (p > 0) # betas
  242        mh = mlag({h}, seq(1, p), e0)
  243            me ~= mh
  244        endif
  245        
  246        dr = rows(me) - rows(vX)
  247        if (dr>0) 
  248        deriv_h = ( zeros(dr, cols(vX)) | vX ) ~ me
  249        else
  250        deriv_h = vX ~ me
  251        endif
  252 
  253        if nmX > 0 
  254        series tmpser = (e>0) ? -1 : 1
  255        matrix sgn = {tmpser}
  256 
  257        matrix mfocs = meanc({e} .* mReg)
  258        matrix dmu = zeros(rows(deriv_h), nmX)
  259 
  260        loop for i = 1 .. q
  261            matrix focs = eeff[,i] .* mlag(mReg, i)
  262            matrix tmpmat = (mlag(sgn, i) + gvec[i]) .* focs
  263            dmu += avec[i] .* tmpmat
  264        endloop
  265        dmu = dmu .* delta
  266 
  267        dmu[1,] = -2*mfocs*sumc(avec|bvec)
  268        deriv_h = dmu ~ deriv_h
  269        endif
  270 
  271        if p > 0
  272        loop for i = 1 .. cols(deriv_h)
  273            series tmpser = deriv_h[,i]
  274            tmpser = filter(tmpser, null, bvec)
  275            deriv_h[,i] = tmpser
  276        endloop
  277        endif       
  278 
  279        if (delta != 2)
  280        deriv_h =  h^(delta/2) .* deriv_h .* (2/delta)
  281        endif
  282    endif
  283    # -------------------------------------------------------------------
  284 
  285    return err
  286 
  287 end function
  288 
  289 /*
  290  * Computes conditional variance and errors for the EGARCH model.
  291  * Returns an error code if any of the parameters is out of its domain. 
  292  */
  293 function scalar egarchFilter(series *y, series *h, series *u, \
  294                  matrix *mReg, matrix *vReg, \
  295                  scalar p, scalar q, matrix *params)
  296 
  297     scalar n_mX = cols(mReg)
  298     scalar n_vX = cols(vReg)
  299     scalar err = 0
  300 
  301     if n_mX == 0
  302     series e = y
  303     else
  304         series e = y - (mReg * params[1:n_mX])
  305     endif
  306 
  307     series u = misszero(e) # Reassigns residuals for the next computation
  308     scalar omegaini = n_mX + 1
  309     scalar base = n_mX + n_vX
  310     if n_vX == 1
  311     series omega = params[omegaini]
  312     else
  313     series omega = vReg * params[omegaini:base]
  314     endif
  315     
  316     series logh = log(var(e))
  317     
  318     # ERRORS
  319     series d = (e>0)
  320     series ae = abs(e)
  321     
  322     if (p == 0) && (q<3)
  323         if q == 1
  324             scalar g = params[base+1]
  325             scalar a = params[base+2]
  326             series tmp = d(-1) ? (g+a) : (g-a)
  327             series logh = omega + tmp*ae(-1)/exp(logh(-1)*0.5)
  328         elif q == 2
  329             scalar g1 = params[base+1]
  330             scalar g2 = params[base+2]
  331             scalar a1 = params[base+3]
  332             scalar a2 = params[base+4]
  333             series tmp1 = d(-1) ? (g1+a1) : (g1-a1)
  334             series tmp2 = d(-2) ? (g2+a2) : (g2-a2)
  335             series logh = omega + tmp1*ae(-1)/exp(logh(-1)*0.5) \
  336             + tmp2*ae(-2)/exp(logh(-2)*0.5)
  337         endif
  338     elif (p == 1) && (q<3)
  339         if q == 1
  340             scalar g = params[base+1]
  341             scalar a = params[base+2]
  342             scalar b = params[base+3]
  343         if (b<1)
  344             series tmp = d(-1) ? (g+a) : (g-a)
  345             series logh = omega + tmp*ae(-1)/exp(logh(-1)*0.5) + b*logh(-1)
  346         else
  347         err = 1
  348         endif
  349 
  350         elif q == 2
  351             scalar g1 = params[base+1]
  352             scalar g2 = params[base+2]
  353             scalar a1 = params[base+3]
  354             scalar a2 = params[base+4]
  355             scalar b  = params[base+5]
  356             series tmp1 = d(-1) ? (g1+a1) : (g1-a1)
  357             series tmp2 = d(-2) ? (g2+a2) : (g2-a2)
  358             series logh = omega + tmp1*ae(-1)/exp(logh(-1)*0.5)) \
  359             + tmp2*ae(-2)/exp(logh(-2)*0.5)) \
  360             + b*logh(-1)
  361         endif
  362     else
  363         string evalstr = "omega"
  364     
  365         loop for i = 1 .. q
  366             scalar a$i = params[base+i]
  367             scalar g$i = params[base+i+p+q]
  368             evalstr += " + a$i*abs(e(-$i)/exp(logh(-$i)*0.5)) + g$i*e(-$i)/exp(logh(-$i)*0.5)"
  369         endloop
  370     
  371         loop for i = 1 .. p
  372             scalar b$i = params[base+p+i]
  373             evalstr += " + b$i*logh(-$i)"
  374         endloop
  375     
  376         series logh = @evalstr
  377     endif
  378 
  379     if (!err) 
  380     if max(logh) > 100
  381         #Overflow check
  382         logh = NA
  383     else
  384         h = exp(logh)
  385     endif
  386     endif
  387 
  388     return err
  389     
  390 end function
  391 
  392 function scalar gfilter(scalar type, 
  393             series *depVar, series *h, series *e, 
  394             matrix *mReg, matrix *vReg, 
  395             scalar p, scalar q, matrix *parameters, 
  396             matrix *DH[null], matrix *DE[null])
  397 
  398     ascore = exists(DH) && exists(DE)
  399 
  400     if type < 7 # aparch
  401     if ascore
  402         err = aparchFilter(&depVar, &h, &e, &mReg, &vReg, p, q, 
  403                    &parameters, has_asymm_fx(type), &DH, &DE)
  404     else
  405         err = aparchFilter(&depVar, &h, &e, &mReg, &vReg, p, q, 
  406                    &parameters, has_asymm_fx(type))
  407     endif
  408 
  409     elif type == 7
  410     err = egarchFilter(&depVar, &h, &e, &mReg, &vReg, p, q, &parameters)
  411     else
  412     err = 1
  413     endif
  414 
  415     return err
  416 
  417 end function
  418 
  419 function matrix do_score(series *de, series *dh, matrix *DE, matrix *DH,
  420              matrix *dd)
  421     ret = {}
  422 
  423     matrix mde = misszero(de)
  424     matrix mdh = misszero(dh)
  425 
  426 #    printf "rows(mde) = %d, rows(mdh) = %d\n", rows(mde), rows(mdh)
  427 #    printf "rows(DE)  = %d, rows(DH)  = %d\n", rows(DE), rows(DH)
  428 
  429 #    printf "%16.9f\n", mde[1:10] ~ DE[1:10,]
  430 #    printf "%16.9f\n", mdh[1:10] ~ DH[1:10,]
  431  
  432     ret = (mde .* DE + mdh .* DH) ~ dd
  433     
  434     return ret
  435 end function
  436 
  437 function scalar do_mle(bundle *model, bool verbose)
  438 
  439     series depVar = model.y
  440     scalar type = model.type
  441 
  442     scalar cdist = model.cdist
  443 
  444     list mlX = model.mlistX
  445     list vlX = model.vlistX
  446 
  447     scalar p = model.p
  448     scalar q = model.q
  449 
  450     scalar nmX  = model.mk
  451     scalar nvX  = model.vk
  452 
  453     scalar err  = 0
  454     series h, e, ll
  455 
  456     mleString = "-"
  457 
  458     if verbose
  459     string mleString += "v"
  460     else
  461     string mleString += "q"
  462     endif
  463     
  464     if model.vcvtype == 0 # robust
  465     string mleString += "r"
  466     elif model.vcvtype == 1 # hessian
  467     string mleString += "h"
  468     # otherwise opg 
  469     endif
  470 
  471     if mleString == "-"
  472     mleString = ""
  473     endif
  474 
  475     fulltheta = model.coeff
  476     inipar = fulltheta
  477     sel = model.active
  478     theta = fulltheta[sel] 
  479 
  480     filtpar_end = rows(fulltheta) - n_cdist_par(cdist)
  481 
  482     mX = model.mX
  483     if model.vX_QR == 1
  484     vX = model.QvX
  485     ini = model.mk+1
  486     fin = ini+model.vk-1
  487     theta[ini:fin] = model.vX_R * theta[ini:fin]
  488     else
  489     vX = model.vX
  490     endif
  491 
  492     set warnings off
  493 
  494     # experimental
  495     set bfgs_toler 1.0e-13
  496 
  497     matrix DH = {}
  498     matrix DE = {}
  499     series dh = NA
  500     series de = NA
  501     matrix dd = {}
  502     matrix score = NA
  503 
  504     if !ascore_ok(type, cdist)
  505     catch mle loglik = ll 
  506             # put the newly estimated values into the full params vector
  507         fulltheta[sel] = theta
  508         filtpar = fulltheta[1:filtpar_end]
  509         matrix dpar = distpar(cdist, fulltheta)
  510         
  511         # FILTER
  512         err = gfilter(type, &depVar, &h, &e, &mX, &vX, p, q, &filtpar)
  513         ll = err ? NA : gig_loglik(&e, &h, cdist, dpar, &de, &dh, &dd)
  514 
  515         params theta
  516     end mle @mleString
  517     err = $error
  518     else
  519     catch mle loglik = ll 
  520             # put the newly estimated values into the full params vector
  521         fulltheta[sel] = theta
  522         filtpar = fulltheta[1:filtpar_end]
  523         matrix dpar = distpar(cdist, fulltheta)
  524         
  525         # FILTER
  526         err = gfilter(type, &depVar, &h, &e, &mX, &vX, p, q, &filtpar, &DH, &DE)
  527         ll = err ? NA : gig_loglik(&e, &h, cdist, dpar, &de, &dh, &dd)
  528         
  529         # SCORE
  530         score = do_score(&de, &dh, &DE, &DH, &dd)
  531         deriv theta = score
  532     end mle @mleString --no-gradient-check
  533     err = $error
  534     endif
  535 
  536   # err = gfilter(type, &depVar, &h, &e, &mX, &vX, p, q, &filtpar, &DH, &DE)
  537   # score = do_score(&de, &dh, &DE, &DH, &dd)
  538   # printf "Score:\n%20.10f\n", sumc(score)
  539     
  540     if (err==0)
  541     matrix crit = {$lnl; $aic; $bic; $hqc}
  542     matrix V = $vcv
  543     else
  544     matrix crit = zeros(4,1)
  545     npar = rows(theta)
  546     matrix V = zeros(npar, npar)
  547     endif
  548 
  549     gig_packResults(&model, err, theta, &h, &e, inipar, V, crit)
  550 
  551     return err
  552 end function