"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020b/addons/SVAR/SVAR_estim_backend.inp" (16 Mar 2020, 8637 Bytes) of package /linux/misc/gretl-2020b.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_estim_backend.inp": 2020a_vs_2020b.

    1 # new file SVAR_estim_backend.inp for SVAR 1.4 (Feb 2019)
    2 # to re-arrange some of the existing functions
    3 # (check LF here)
    4 
    5 ###########################
    6 
    7 function scalar base_est(bundle *SVARobj)
    8 
    9 /*
   10    takes a SVAR object and adds the basic VAR estimates
   11    to the bundle; returns an errcode (virginity check).
   12     */
   13 
   14     scalar step = SVARobj.step
   15     err = (step>0)
   16 
   17     if err
   18         printf "Base estimation done already!\n"
   19     else
   20         matrix mY  = SVARobj.Y
   21         scalar p   = SVARobj.p
   22         scalar n   = SVARobj.n
   23         scalar k   = SVARobj.k
   24 
   25         scalar T   = SVARobj.T
   26 
   27         matrix E
   28         matrix mreg = SVARobj.X ~ mlag(mY, seq(1,p))
   29         matrix olspar = mols(mY[p+1:T,], mreg[p+1:T,], &E)
   30         matrix Sigma = (E'E) ./ (T - (n*p + k)) # was:  / df
   31 
   32         matrix SVARobj.VARpar = transp(olspar[k+1:k+n*p,])
   33         matrix SVARobj.mu = (k>0) ? olspar[1:k,] : {} # was: d
   34         matrix SVARobj.Sigma = Sigma
   35         matrix SVARobj.E = E
   36         SVARobj.step = 1
   37         SVARobj.LL0 = T * (n* log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
   38     endif
   39 
   40     return err
   41 
   42 end function
   43 
   44 ####################
   45 
   46 function scalar vecm_est(bundle *SVARobj)
   47 
   48 /*
   49    We can't afford to be too flexible here, and the intended usage
   50    is as follows.
   51 
   52    We assume the model has already been declared as a SVEC (type=4)
   53    model. We also assume that the cointegration properties (beta mat plus
   54    deterministics) have already been set up via the SVAR_coint() function,
   55    so that we already have beta and alpha available as "jbeta" and
   56    "jalpha", respectively. Finally, we take care of proper treatment of
   57    deterministics, via the scalar "jcase" (1 to 5).
   58     */
   59 
   60     # --- preliminary checks
   61 
   62     err = (SVARobj.step > 0)
   63 
   64     if err
   65         printf "Base estimation done already!\n"
   66         return err
   67     endif
   68 
   69     err = (SVARobj.type != 4)
   70 
   71     if err
   72         printf "Wrong model type!\n"
   73         return err
   74     endif
   75 
   76     # --- grab stuff from the bundle
   77 
   78     matrix mY     = SVARobj.Y
   79     matrix jbeta   = SVARobj.jbeta
   80     matrix jalpha  = SVARobj.jalpha
   81     scalar p      = SVARobj.p
   82     scalar n      = SVARobj.n
   83     scalar k      = SVARobj.k
   84     scalar r      = cols(jbeta)
   85     scalar dcase   = SVARobj.jcase
   86     scalar ols_al = rows(jalpha) == 0
   87 
   88     scalar T      = SVARobj.T
   89 
   90     # --- first thing: do we have a pre-set value for alpha?
   91 
   92     matrix dY = diff(mY)
   93     matrix dep = dY[p+1:T,]
   94     matrix E = {}
   95     ng = n * (p-1) # number coming from Gammas /lagged diffs (p>0 for SVEC)
   96 
   97     # deterministics
   98     matrix mreg = vecm_det(T, dcase)
   99 
  100     # ECM terms
  101     matrix ECM = mlag(mY * jbeta[1:n,], 1)
  102     if dcase == 2
  103         ECM = ECM .+ jbeta[n+1, ]
  104     elif dcase == 4
  105         ECM += seq(1,T)'jbeta[n+1, ]
  106     endif
  107 
  108     if ols_al
  109         # alpha must be estimated together with the rest of the VECM
  110         matrix mreg ~= SVARobj.X ~ ECM
  111     else
  112         matrix dep -= (ECM[p+1:T,] * jalpha')
  113         matrix mreg ~= SVARobj.X
  114     endif
  115 
  116     # extra lags
  117     if p > 1
  118         mreg ~= mlag(dY, seq(1,p-1))
  119     endif
  120 
  121     # trim the first rows
  122     if rows(mreg)>0
  123         mreg = mreg[p+1:T,]
  124         matrix olspar = mols(dep, mreg, &E)
  125     else
  126         matrix olspar = {}
  127         E = dep
  128     endif
  129 
  130     df = T - (n*p + k - (n-r)) # FIXME: aren't the unrestr. determ. terms missing?
  131     # (but apparently df is never used again...)
  132     matrix Sigma = (E'E)./T
  133 
  134     # --- construct the various matrices required later
  135 
  136     rp = rows(olspar)
  137 
  138     # alpha first (should be before the gammas, if estimated)
  139     if ols_al
  140         jalpha = olspar[rp-ng-r+1 : rp-ng,]'
  141     endif
  142 
  143     # exogenous
  144 
  145     if dcase == 1
  146         matrix mu = {}
  147         scalar nd = 0
  148     elif dcase == 2
  149         matrix mu = jbeta[n+1,] * jalpha'
  150         scalar nd = 0
  151     elif dcase == 3
  152         matrix mu = olspar[1,]
  153         scalar nd = 1
  154     elif dcase == 4
  155         matrix mu =  olspar[1,] |(jbeta[n+1,] * jalpha')
  156         scalar nd = 1
  157     elif dcase == 5
  158         matrix mu = olspar[1:2,]
  159         scalar nd = 2
  160     endif
  161 
  162     if k > 0
  163         mu = mu | olspar[nd + 1 : nd + k,] # was: mu ~ olspar[sel,]
  164     endif
  165 
  166     /*
  167        companion form in levels
  168     */
  169     Pi = jalpha * jbeta[1:n,]'
  170 
  171     if p > 1
  172         # the Gammas are always at the back:
  173         ini = rp - ng + 1
  174 
  175         matrix A = olspar[ini: ,] | zeros(n,n)  # was: ini:fin
  176         A += I(n) + Pi' | -olspar[ini: ,]
  177     else
  178         matrix A = I(n) + Pi'
  179     endif
  180 
  181     matrix SVARobj.VARpar = A'
  182     matrix SVARobj.mu = mu
  183     matrix SVARobj.Sigma = Sigma
  184     matrix SVARobj.E = E
  185     matrix SVARobj.jalpha = jalpha
  186     SVARobj.step = 1
  187     SVARobj.LL0 = T*(n*log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
  188 
  189     return err
  190 
  191 end function
  192 
  193 ##########
  194 
  195 function scalar VARloglik(scalar T, matrix Sigma, matrix *C[null])
  196     # computes the (concentrated) loglikelihood for a VAR
  197     # the matrix C (such that in theory CC' = Sigma) may be null,
  198     # for just-identified models
  199 
  200     n = rows(Sigma)
  201     ll = n * 1.83787706640934548355 # ln(2*$pi)
  202 
  203     if !exists(C) # just-identified model
  204         ll += log(det(Sigma)) + n
  205     else
  206         matrix KK = invpd(C*C')
  207         ll += -log(det(KK)) + tr(KK*Sigma)
  208     endif
  209 
  210     return -(T/2) * ll
  211 end function
  212 
  213 ##########
  214 
  215 function scalar loglik(matrix *theta, matrix *dat, scalar modeltype)
  216     # modeltype < 0->C; modeltype>=0: AB (contains free elements of a)
  217     matrix Sigma = {}
  218     matrix Ss = {}
  219     unscramble_dat(&dat, &Sigma, &Ss)
  220     if modeltype < 0
  221         matrix C = mat_exp(theta, Ss, 0)
  222     else
  223         p1 = modeltype
  224         p2 = rows(theta)
  225         matrix aSs = Ss[,1:p1+1]
  226         matrix bSs = Ss[,p1+2:p2+2]
  227         matrix A B
  228         ABmat_exp(theta, aSs, bSs, &A, &B)
  229         # was: matrix C = B/A
  230         matrix C = A\B
  231     endif
  232 
  233     matrix KK = invpd(C*C')
  234     ll = det(KK) # should always be positive
  235     ll = (ll<=0) ? NA : -0.5 * (tr(KK*Sigma) - log(ll))
  236 
  237     return ll
  238 end function
  239 
  240 ##############
  241 
  242 
  243 function matrix InfoMat(matrix CorB, matrix S, matrix *A[null])
  244 /*
  245    merged from InfoMatC and InfoMatAB
  246    First case C model (A is null), latter AB model.
  247     */
  248     matrix tmp = !exists(A) ? I(rows(CorB)) : (A\CorB) | -I(rows(A))
  249     tmp = S' (tmp ** inv(CorB'))
  250     return tmp * N2_ify(tmp')
  251 end function
  252 
  253 
  254 function matrix coeffVCV(matrix S, matrix *rhs, matrix *lhs[null])
  255     # C or AB model
  256     matrix IM = !exists(lhs) ? InfoMat(rhs, S) : InfoMat(rhs, S, &lhs)
  257     # (should be correct with new InfoMat)
  258 
  259 
  260     # quick-dirty check for singularity
  261     if rcond(IM)>1.0e-10
  262         matrix iIM = invpd(IM)/$nobs
  263     else
  264         matrix evec
  265         l = eigensym(IM, &evec)
  266         printf "\n\nInformation matrix is not pd!!\n\n%12.5f\n", IM
  267         printf "Eigenvalues:\n%12.5f\n", l
  268         matrix e = (l .< 1.0e-07)
  269         printf "Troublesome eigenvectors:\n%12.5f\n", selifc(evec, e')
  270         printf "S:\n%4.0f\n", S
  271         matrix iIM = zeros(rows(IM), cols(IM))
  272     endif
  273 
  274     return qform(S,iIM)
  275 end function
  276 
  277 ################
  278 
  279 function void SVAR_est_printout(bundle *mod)
  280     # scalar type = mod.type
  281 
  282     printf "Optimization method = "
  283     strings os = defarray("BFGS (numerical)", "BFGS (analytical)", \
  284       "Newton-Raphson (numerical)", "Newton-Raphson (analytical score)", \
  285       "Scoring algorithm")
  286     printf "%s\n", os[mod.optmeth + 1]  # Because optmeth starts at 0 here
  287 
  288     printf "Unconstrained Sigma:\n%12.5f\n", mod.Sigma
  289 
  290     if mod.type < 3 || mod.type == 4 # plain or C-model, or (Jan 2018) SVEC
  291         # Standard C matrix
  292         printStrMat(mod.C, mod.vcv, "C")    # was mod.S1
  293 
  294         # Long-run matrix
  295         if rows(mod.Rd1l) || mod.calc_lr || mod.type == 4
  296             # either long-run restr., or user wish, or SVEC
  297             # (this could be much embellished, probably)
  298             printf "Estimated long-run matrix"
  299             if rows(mod.Rd1l)
  300                 printf " (restricted)"
  301             endif
  302             printf "\n"
  303 
  304             matrix longrun = mod.lrmat
  305             # round to zero for printout (also do it in general?)
  306             longrun = (abs(longrun) .< 1e-15) ? 0.0 : longrun
  307             print longrun
  308         endif
  309 
  310     elif mod.type == 3  # AB, new <Sven>
  311         n2 = (mod.n)^2
  312         if mod.ka > 0
  313             printStrMat(mod.S1, mod.vcv[1:n2, 1:n2], "A")
  314         endif
  315 
  316         if mod.kb > 0
  317             printStrMat(mod.S2, mod.vcv[n2+1 : 2*n2, n2+1 : 2*n2], "B")
  318         endif
  319     endif
  320 
  321     printf "  Log-likelihood = %g\n", mod.LL1
  322     if inbundle(mod, "LRoid")   # over-id test was done
  323         printf "  Overidentification LR test = %g (%d df, pval = %g)\n\n", \
  324           mod.LRoid[1], mod.LRoid[2], mod.LRoid[3]
  325     endif
  326 
  327 end function