"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/addons/SVAR/SVAR_Cfuncs.inp" (15 Sep 2020, 5848 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 "SVAR_Cfuncs.inp": 2020d_vs_2020e.

    1 function void scoreC(matrix *ret, matrix *theta, matrix *dat)
    2     matrix Sigma Ss
    3     unscramble_dat(&dat, &Sigma, &Ss)
    4     scalar n = rows(Sigma)
    5     scalar npar = cols(Ss) - 1
    6 
    7     matrix C = mat_exp(theta, Ss) # was: , 0)
    8     matrix S = Ss[,1:npar]
    9 
   10     matrix iC = inv(C)
   11     matrix ret = iC' (qform(iC,Sigma) - I(n))
   12 
   13     ret = vec(ret)'S
   14 end function
   15 
   16 function matrix C1mat( const matrix A, bool VECM[0],
   17                       const matrix jalpha[null], const matrix jbeta[null])
   18     # (Jan 2018: used to be pointer args, coming from a time where it was necessary
   19     # for optional args; now changed to non-pointers)
   20 
   21     /*
   22        computes C(1) out of the autoregressive matrices; jalpha and jbeta
   23        are alpha and beta, which are used only if VECM is nonzero, that
   24        is under cointegration
   25 
   26        (Remark Sven: here jbeta was assumed to be (n,r)
   27        without coeffs for restr. terms
   28        Different from jbeta (n+1,r) in other contexts (?).
   29     We now (Jan 2018) make sure only rows 1:n are used. */
   30 
   31     scalar n = rows(A)
   32     scalar p = cols(A) / n
   33     matrix tmp = mshape(vec(A), n*n, p)
   34 
   35     if VECM == 0
   36         tmp = mshape(sumr(tmp), n, n)
   37         matrix ret = inv(I(n) - tmp)
   38     else # cointegrated
   39         if !exists(jalpha) || !exists(jbeta)
   40             funcerr "Need cointegration params for C1 in VECM!"
   41         endif
   42 
   43         matrix aperp = nullspace(jalpha')
   44         matrix bperp = nullspace(jbeta[1:n, ]')
   45         tmp = mshape(tmp*seq(1,p)', n, n)
   46         matrix ret = bperp * inv(aperp'tmp*bperp) * aperp'
   47     endif
   48 
   49     return ret
   50 end function
   51 
   52 function matrix init_C(matrix Sigma, matrix Rd)
   53     matrix Ss = imp2exp(Rd)
   54     scalar k = cols(Ss)
   55     if k == 1
   56         # nothing to estimate
   57         ret = {}
   58     else
   59         scalar n = rows(Sigma)
   60 
   61         matrix S = (k>1) ? Ss[,1:k-1] : {}
   62         matrix s = Ss[,k]
   63         matrix K = cholesky(Sigma)
   64         matrix bigmat = (K ** I(n))
   65 
   66         matrix ret = mols(vec(Sigma) - bigmat*s, bigmat*S)
   67     endif
   68 
   69     return ret
   70 end function
   71 
   72 function void PseudoHessC(matrix *H, matrix *theta, matrix *dat)
   73     matrix Sigma Ss
   74     unscramble_dat(&dat, &Sigma, &Ss)
   75     scalar n = rows(Sigma)
   76     scalar npar = cols(Ss) - 1
   77     matrix S = Ss[,1:npar]
   78 
   79     # printf "PseudoHessC\n%10.4f\n%4.1f\n", theta, Ss
   80 
   81     matrix C = mat_exp(theta, Ss) # was: , 0)
   82     H = InfoMat(C, S) # was InfoMatC
   83 end function
   84 
   85 function matrix estC(matrix *theta, matrix Sigma, matrix Rd, matrix *vcv[null],\
   86                      scalar *errcode, int method, scalar verbose[1])
   87 
   88     matrix Ss = imp2exp(Rd)
   89     scalar n = rows(Sigma)
   90     if cols(Ss) == 1
   91         # nothing to estimate
   92         matrix C = mshape(Ss, n, n)
   93         if exists(vcv)
   94             vcv = zeros(n*n,n*n)
   95         endif
   96         errcode = 0
   97         return C
   98     endif
   99 
  100     scalar SCALE = 0 # EXPERIMENTAL
  101     scalar npar = rows(theta)
  102 
  103     if SCALE == 1
  104         printf "Scale!\n"
  105         matrix s = sqrt(diag(Sigma))
  106         matrix sSig = Sigma ./ (s*s')
  107     else
  108         matrix sSig = Sigma
  109     endif
  110 
  111     # printf "estC\n%10.4f\n%4.1f\n", theta, Ss
  112 
  113     matrix tmp = mat_exp(theta, Ss) # was: , 1)
  114 
  115     if verbose > 2
  116         /* obsolete ? */
  117         printf "check within estC -- before estimation\n"
  118         check_const(tmp, Rd)
  119     endif
  120     matrix dat = vec(sSig) ~ Ss
  121 
  122     matrix g H
  123     if verbose > 1
  124         set max_verbose 1
  125     else
  126         set max_verbose 0
  127     endif
  128     err = 1
  129     iters = 0
  130     #set bfgs_toler 1.0e-03
  131 
  132     matrix theta0 = theta
  133     errcode = 0
  134 
  135     loop while (err==1 && iters<100) --quiet
  136         if method == 0
  137             catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, -1)")
  138             errcode = $error
  139             scoreC(&g, &theta, &dat)
  140         elif method == 1
  141             catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, -1)",
  142               "scoreC(&g, &theta, &dat)")
  143             errcode = $error
  144         elif method == 2
  145             catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)")
  146             errcode = $error
  147             scoreC(&g, &theta, &dat)
  148         elif method == 3
  149             catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)",
  150               "scoreC(&g, &theta, &dat)")
  151             errcode = $error
  152         elif method == 4
  153             catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)",
  154               "scoreC(&g, &theta, &dat)",
  155               "PseudoHessC(&H, &theta, &dat)")
  156             errcode = $error
  157         endif
  158 
  159         if errcode>0
  160             printf "errcode = %d\n", errcode
  161         endif
  162         scalar crit = maxr(abs(g))
  163         err = (crit > 1.0e-4)
  164 
  165         if err==1
  166             iters++
  167             theta = 0.1*mnormal(npar,1) + theta0
  168             if verbose>1
  169                 printf "Iter %3d: Restarting... ll = %12.7f, crit = %16.10f\n", \
  170                   iters, ll, crit
  171                 printf "theta = %10.5f grad = %10.5f\n", theta', g
  172             endif
  173         endif
  174 
  175     endloop
  176 
  177     scalar crit = maxr(abs(g))
  178     warn = (crit > 1.0e-1)
  179 
  180     if !err
  181         if (iters > 0) && (verbose > 1)
  182             printf "Converged after %d restarts\n", iters
  183         endif
  184 
  185         # printf "estC(2)\n%10.4f\n%4.1f\n", theta, Ss
  186         matrix C = mat_exp(theta, Ss) # was: , 1)
  187 
  188         if SCALE == 1
  189             C = C .* s'
  190         endif
  191         if verbose > 1
  192             printf "Estimated C matrix:\n%12.5f", C
  193         endif
  194 
  195         if exists(vcv)
  196             vcv = coeffVCV(Ss[,1:npar], &C)
  197         endif
  198 
  199         if verbose > 1
  200             matrix Info = InfoMat(C, Ss[,1:npar])   # was InfoMatC
  201             printf "estC : Info = \n%14.10f\n", Info
  202             printf "estC : score = \n%14.10f\n", g
  203         endif
  204 
  205     else
  206         if verbose > 1
  207             printf "No convergence! :-( \n%12.6f\n", theta' | g
  208         endif
  209         matrix C = {}
  210         errcode = 1
  211     endif
  212 
  213     return C
  214 end function