"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/addons/SVAR/SVAR_IRF.inp" (11 Apr 2020, 8787 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_IRF.inp": 2020a_vs_2020b.

    1 function void doIRF(bundle *SVARobj)
    2 /*
    3    constructs the structural VMA representation. Note
    4    that the companion matrix is never used explicitly;
    5 
    6    The output is not returned by the function, but rather
    7    put into the bundle under the "IRFs" key.
    8     */
    9     scalar type = SVARobj.type
   10     matrix varA = SVARobj.VARpar
   11     scalar H = SVARobj.horizon + 1
   12     scalar n = SVARobj.n
   13 
   14     if (type == 1) || (type == 2) || (type == 4)
   15         matrix C = SVARobj.C    # was: SVARobj.S1
   16     elif type == 3
   17         if inbundle(SVARobj, "C")   # maybe not yet computed
   18             matrix C = SVARobj.C
   19         else
   20             matrix C = SVARobj.S1 \ SVARobj.S2
   21         endif
   22     endif
   23 
   24     matrix ret = zeros(H,n*n)
   25     scalar np = SVARobj.p * n
   26     matrix tmp = I(np)
   27     matrix prd = zeros(np,np)
   28 
   29     loop i=1..H --quiet
   30         ret[i,] = vec(tmp[1:n,1:n] * C)'
   31         if (np>n)
   32             prd[n+1:np, ] = tmp[1:np-n, ]
   33         endif
   34         prd[1:n,] = varA * tmp
   35         tmp = prd
   36     endloop
   37 
   38     if SVARobj.ncumul > 0
   39         # The following code is now done in SVAR_cumulate
   40         # once and for all:
   41         # matrix to_cum = SVARobj.cumul
   42         # tmp = zeros(n,n)
   43         # tmp[to_cum,] = 1
   44         # sel = selifr(transp(seq(1,n*n)), vec(tmp))
   45         ret[, SVARobj.cumsel] = cum(ret[, SVARobj.cumsel]) # .cumsel was sel
   46     endif
   47 
   48     matrix SVARobj.IRFs = ret
   49 end function
   50 
   51 ##############
   52 
   53 function matrix FEVD(bundle *SVARobj, int drawix[0::0])
   54     # (drawix only meant for the set id case (type 10))
   55 
   56     n = SVARobj.n
   57     h = SVARobj.horizon + 1
   58 
   59     if SVARobj.type == 10 
   60         # In the set id case in each accepted draw the impulse
   61         # responses are already stored as "irfs"; however, the format 
   62         # there is an array of matrices. 
   63 
   64         errchkSRhisto(&SVARobj, drawix)
   65         
   66         # allow drawix to override the setting in the bundle
   67         whichdraw = drawix ? drawix : SVARobj.bestdraw
   68         bundle pickdraw = SVARobj.acc_draws[whichdraw]
   69 
   70         if cols(pickdraw.irfs[1]) != n
   71             funcerr "partial id not supported for FEVD"
   72         elif h != nelem(pickdraw.irfs)
   73             funcerr "horizon mismatch"
   74         endif
   75 
   76         if !inbundle(pickdraw, "IRFs") # maybe have already been added there
   77             putIrf_to_accdraw(&SVARobj, whichdraw) 
   78             
   79             # matrix IRFs = zeros(h, n*n)
   80             # loop ix = 1..h -q
   81             #     IRFs[ix, ] = vec(pickdraw.irfs[ix])'
   82             # endloop
   83             ## copy to origin
   84             # matrix SVARobj.acc_draws[whichdraw].IRFs = IRFs
   85         # else 
   86         #     matrix IRFs = pickdraw.IRFs
   87         endif
   88         matrix IRFs = SVARobj.acc_draws[whichdraw].IRFs
   89 
   90     else # standard non-SR model 
   91         if drawix > 0
   92             print "Warning: 'drawix' arg meaningless for standard SVAR, ignoring"
   93         endif
   94         matrix IRFs = SVARobj.IRFs
   95     endif
   96 
   97     matrix ret = zeros(h, n*n)
   98     ctmp = cum(IRFs .* IRFs)
   99 
  100     loop i = 1..h --quiet
  101         tmp = mshape(ctmp[i,],n,n)'
  102         ret[i,] = vec(tmp ./ sumc(tmp))'
  103     endloop
  104 
  105     return ret
  106 end function
  107 
  108 
  109 ### The functions GetShock() and SVAR_getshock aren't called anywhere, 
  110 #### but they're public and meant to be used by the user.
  111 ## GetShock may be deprecated in the future in favor of SVAR_getshock.
  112 
  113 function series SVAR_getshock(bundle *mod, string sname[null], 
  114                               int drawix[0::0])
  115     # This is a wrapper to provide a nicer interface, using the 
  116     # name of the shock instead of the number.
  117     # Default (as in GetShock) is to use the first shock.
  118 
  119     s_ix = !exists(sname) ? 1 : strpos_allin(mod.snames, sname)
  120 
  121     return GetShock(&mod, s_ix, drawix)
  122 end function 
  123 
  124 function series GetShock(bundle *SVARobj, int i[1::1], int drawix[0::0])
  125     /*
  126     Produces the series corresponding to the historical shock 
  127     realizations associated with the point estimates of the model
  128     (and IRFs).
  129     For set identification (sign restrictions) there is no point
  130     estimate; however, we support that 
  131     the user picks one of the accepted draws and then the shock series 
  132     is based on that particular model draw.
  133     # (drawix only meant for the set id case (type 10))
  134     */
  135 
  136     series ret = NA
  137     type = SVARobj.type
  138     matrix B10 = {} # to be filled in type 10
  139 
  140     ## some error checks ##
  141     if type > 4 && type != 10
  142         printf "Given type %d\n", type
  143         funcerr "Unknown model type"
  144     
  145     elif i > SVARobj.n
  146         printf "Chosen shock index: %d\n", i
  147         funcerr "Shock index out of range"
  148         
  149     elif type != 10 && drawix > 0
  150         print "Warning: 'drawix' arg meaningless for standard SVAR"
  151     elif type == 10 
  152         errchkSRhisto(&SVARobj, drawix)
  153     endif
  154 
  155     ## get the C matrix (and then the inv) ##
  156     if (type == 1) || (type == 2) || (type == 4)
  157         matrix C = SVARobj.C 
  158 
  159     elif type == 3
  160         if inbundle(SVARobj, "C")   # maybe not yet computed
  161             matrix C = SVARobj.C
  162         else
  163             matrix C = SVARobj.S1 \ SVARobj.S2
  164         endif
  165 
  166     elif type == 10 # set id 
  167         # allow drawix to override the setting in the bundle
  168         whichdraw = drawix ? drawix : SVARobj.bestdraw
  169         bundle pickdraw = SVARobj.acc_draws[whichdraw]
  170         matrix C = pickdraw.irfs[1] # impact effect is C
  171         B10 = pickdraw.B
  172 
  173         if cols(C) < SVARobj.n
  174             funcerr "partial id not supported for shock retrieval"
  175         endif
  176     endif
  177 
  178     matrix iC = inv(C')
  179     
  180     matrix resids = muVARparE_mayberedr(SVARobj, B10)[3]
  181 
  182     ## construct the wanted series ##
  183     extra = $nobs - rows(resids)
  184     matrix tmp = {}
  185     if extra > 0
  186         set warnings off
  187         tmp = ones(extra,1) .* NA
  188     endif
  189 
  190     tmp |= resids * iC[,i]
  191     ret = tmp
  192 
  193     snames = SVARobj.snames # strings array?
  194     string vlab = snames[i]
  195 
  196     setinfo ret --description="@vlab"
  197 
  198     return ret
  199 end function
  200 
  201 #######################
  202 
  203 function list SVAR_HD(bundle *mod, string vname[null], 
  204                       int drawix[0::0])
  205     # wrapper around SVAR_hd to use a string interface for the 
  206     # variable 
  207     v_ix = !exists(vname) ? 1 : strpos_allin(mod.Ynames, vname)
  208     return SVAR_hd(&mod, v_ix, drawix)  
  209 end function
  210 
  211 
  212 function list SVAR_hd(bundle *Mod, int nv[1::1], int drawix[0::0])
  213     # historical decomposition
  214     # (drawix only meant for the set id case (type 10))
  215 
  216     list ret = null
  217     loop foreach i n p t1 t2 type T k --quiet
  218         scalar $i = Mod.$i
  219     endloop
  220     matrix B10 = {} # to be filled for type 10
  221 
  222     if nv > n
  223         printf "Hm. There are %d variables in the model. ", n 
  224         printf "Chosen shock index: %d\n", nv
  225         funcerr "Shock index out of range"
  226         ## (further range check for SR partial id below) ##
  227     endif
  228     
  229     # Prepare the set id case
  230     if type == 10 
  231         errchkSRhisto(&Mod, drawix)
  232         # allow drawix to override the setting in the bundle
  233         whichdraw = drawix ? drawix : Mod.bestdraw
  234         bundle pickdraw = Mod.acc_draws[whichdraw] 
  235         matrix B10 = pickdraw.B  
  236     endif
  237 
  238     # The following might be redrawn in type10/Bayesian
  239     matrices muVARparE = muVARparE_mayberedr(Mod, B10)
  240 
  241     # compute the exogenous part
  242     if type < 4
  243         matrix m = Mod.X * Mod.mu
  244 
  245     elif type == 4
  246         # here we have to take into account the "5 cases"
  247         dcase = Mod.jcase
  248         # T     = Mod.T
  249         matrix mreg = (dcase == 1) ? {} : ones(T,1)
  250         if dcase > 3
  251             mreg ~= seq(1,T)'
  252         endif
  253 
  254         matrix m = (mreg ~ Mod.X) * Mod.mu
  255 
  256     elif type == 10
  257         matrix m = Mod.X * muVARparE[1]
  258     endif 
  259 
  260     # grab the C matrix
  261     if (type == 1) || (type == 2) || (type == 4)
  262         matrix C = Mod.C
  263 
  264     elif type == 3
  265         if inbundle(Mod, "C")
  266             matrix C = Mod.C
  267         else
  268             matrix C = Mod.S1 \ Mod.S2
  269         endif
  270 
  271     elif type == 10
  272         matrix C = pickdraw.irfs[1] # impact effect is C
  273         if cols(C) < Mod.n
  274             funcerr "partial id not supported for historical decomp"
  275         endif
  276     endif
  277 
  278     matrix iC = inv(C)
  279     strings Ynames = Mod.Ynames
  280     strings snames = Mod.snames
  281     string yn = Ynames[nv]
  282 
  283     smpl t1 t2
  284     if cols(m)>0
  285         Xdet = varsimul(muVARparE[2], m[p+1:,], Mod.Y[1:p,]) # was VARpar
  286     else
  287         Xdet = varsimul(muVARparE[2], zeros(T-p, n), Mod.Y[1:p,]) # was Mod.T
  288     endif
  289 
  290     ret += genseries( sprintf("hd_%s_det", yn), Xdet[,nv])
  291 
  292     # the structural shocks
  293     matrix U = muVARparE[3] * iC' # was E
  294     rotVARpar = iC * muVARparE[2] * (I(p) ** C)
  295 
  296     loop i = 1..n --quiet
  297         a = (seq(1,n) .= i)
  298         W = varsimul(rotVARpar, U .* a, zeros(p,n)) * C'
  299         ret += genseries(sprintf("hd_%s_%s", yn,
  300                  fixname(snames[i])), W[,nv])
  301     endloop
  302 
  303     return ret
  304 end function