Am 02.07.2018 um 19:56 schrieb Sven Schreiber: > Am 02.07.2018 um 19:00 schrieb Andreas Zervas: >> I am attaching a test script. For case 3 it works, but not for case >> 4, and as Sven said, nor for case 2. >> > Ok, thanks. If there are no further complications, I think I should be > able to deal with this tomorrow.
It seems my hunch was correct, it was simply a bug in the matrix multiplication order, so nobody before used an SVEC with a rank higher than one and with the restricted deterministics setup -- or at least nobody before cared to report the problem... In vecm_est() the following two lines needed to be changed: matrix mu = jbeta[n+1,] * jalpha' # was: jalpha * jbeta[n+1,] ... matrix mu = olspar[1,] | (jbeta[n+1,] * jalpha') # was: (jalpha * jbeta[n+1,])' It works for me, but I'm attaching a new gfn for testing -- but beware this is not the official updated release, for example the sample files and the help is missing there. Also attached is a slightly generalized test script. I'm not currently at a machine where the git access is set up, so pushing this to git would have to wait until tonight (unless Jack does it first). cheers, sven
|
Am 02.07.2018 um 19:56 schrieb Sven
Schreiber:
It seems my hunch was correct, it was simply a bug in the matrix multiplication order, so nobody before used an SVEC with a rank higher than one and with the restricted deterministics setup -- or at least nobody before cared to report the problem... In vecm_est() the following two lines needed to be changed: matrix mu = jbeta[n+1,] * jalpha' # was: jalpha * jbeta[n+1,] ... matrix mu = olspar[1,] | (jbeta[n+1,] * jalpha') # was: (jalpha * jbeta[n+1,])' It works for me, but I'm attaching a new gfn for testing -- but beware this is not the official updated release, for example the sample files and the help is missing there. Also attached is a slightly generalized test script. I'm not currently at a machine where the git access is set up, so pushing this to git would have to wait until tonight (unless Jack does it first). cheers, sven |
<?xml version="1.0" encoding="UTF-8"?>
<gretl-functions>
<gretl-function-package name="SVAR" needs-time-series-data="true" minver="2017a" lives-in-subdir="true">
<author email="r.lucche...@univpm.it">Riccardo "Jack" Lucchetti and Sven Schreiber</author>
<version>1.36</version>
<date>2018-07-03</date>
<description>Structural VARs</description>
<tags>C32</tags>
<label>Structural VARs</label>
<help>
pdfdoc:SVAR.pdf
</help>
<data-files count="1">
examples </data-files>
<gretl-function name="SVAR_setup" type="bundle">
<params count="5">
<param name="type_string" type="string"/>
<param name="lY" type="list"/>
<param name="lX" type="list" optional="true"/>
<param name="varorder" type="int" min="1" default="1"/>
<param name="dcase" type="int" min="1" max="5" default="3"/>
</params>
<code>/*
This creates a bundle holding the info on the SVAR; this will be
filled in 3 steps: This function inserts the initial info: sample
size, data and so on. Then, more stuff will have to be added later:
the VAR estimates in packed form (see below), then the SVAR estimates.
The scalar "step" keeps track of the stage we're at.
(maybe not needed any more)
The bootstrap does _not_ live here.
*/
scalar type = modeltype(type_string)
if type == 4 # SVEC
# SVEC treats restricted deterministic terms specially as per Johansen,
# so we'll drop them if present,
# any other exogenous regressors are under the user's responsibility
# (eg centred dummies)
# the following currently doesn't work in gretl w.r.t. the time trend,
# but fails silently (but const works)
list lX -= const time
# therefore check for remaining collinearity
genr time
list lcheck = lX const time
matrix mcheck = { lcheck }
if rank(mcheck'mcheck) < nelem(lcheck)
funcerr "Remove constant and trend terms from X list!"
endif
endif
scalar n = nelem(lY)
scalar k = nelem(lX)
scalar n2 = n*n
bundle ret = null
/*-----------------------------------------------------------------------*/
/* general info */
/*-----------------------------------------------------------------------*/
/* type goes from 1 to 4 (plain, "C", "AB" or "SVEC") */
ret.type = type
/* the step we're at: 0 = no estimation done, 1 = VAR only, 2 = SVAR */
ret.step = 0 # maybe not needed any more
matrix mreg = set_default_dimensions(&ret, lY, lX, varorder)
/* the actual data */
matrix ret.Y = mreg[,1:n]
matrix ret.X = (k>0) ? mreg[, n+1 : n+k] : {}
# variable names
strings ret.Ynames = strsplit(strsub(varname(lY),","," ")) # note: string array
strings ret.Xnames = strsplit(strsub(varname(lX),","," ")) # needed?
/*
The constraint matrices.
"Rd1" contains short-run constraints on B (and therefore C in non-AB models);
"Rd1l" contains long-run constraints on C (not supported in AB models);
"Rd0" contains short-run constraints on A in AB;
note that "Rd1l" and "Rd0" were both "aux" in previous versions.
Initially, they are empty. Except for the "plain" model, it's up to
the user to fill them up later, via SVAR_restrict() or by hand
*/
matrix ret.Rd1 = (type==1) ? cholRd(n) : {}
matrix ret.Rd1l = {}
if type == 3
matrix ret.Rd0 = {}
endif
if type == 4
/* SVEC model: cointegration stuff */
matrix ret.jalpha = {}
matrix ret.jbeta = {}
ret.jcase = 0
endif
/*
Optimisation
*/
ret.optmeth = 4 # default = scoring
ret.simann = 0 # not for now
/*
Information for cumulating/normalizing IRFs
*/
ret.ncumul = 0
matrix ret.cumul = {}
ret.normalize = 0
/*
Bootstrap
*/
ret.nboot = 0
ret.boot_alpha = -1
matrix ret.bootdata = {}
ret.biascorr = 0
/* Names for shocks */
# per default, we borrow variable names, following tradition
strings ret.snames = ret.Ynames
/* Don't calculate the long-run matrix by default.
(This does not apply to the case of long-run restrictions.) */
ret.calc_lr = 0
return ret
</code>
</gretl-function>
<gretl-function name="SVAR_restrict" type="scalar">
<params count="5">
<param name="b" type="bundleref"/>
<param name="code" type="string"/>
<param name="r" type="int"/>
<param name="c" type="int" default="0"/>
<param name="d" type="scalar" default="0"/>
</params>
<code># c gets a default so that it can be omitted with Adiag, Bdiag (?)
# the d default is also natural
type = b.type
n = b.n
# check input for implemented restriction code
if strstr("C lrC A B Adiag Bdiag", code) == "" # code unknown
print "Unknown code in SVAR_restrict."
return 2
# check for input mismatch
elif (code=="C" || code=="lrC") && !(type==1 || type==2 || type==4)
print "C type restriction but not a C model."
return 1
elif (strstr("A B Adiag Bdiag", code) != "") && (type!=3)
print "AB type restriction but not an AB model."
return 1
# check for unsupported case
elif code=="lrC" && type==3
print "Long-run restrictions only supported in C model."
return 3 # (Which code to return here? <Sven>)
endif
# if no input error, proceed with this:
err = 0
if code == "C"
matrix Rd = b.Rd1
err = add_constraint(&Rd, n, r, c, d)
if !err
b.Rd1 = Rd
endif
elif code == "lrC"
matrix Rd = b.Rd1l
err = add_constraint(&Rd, n, r, c, d)
if !err
b.Rd1l = Rd
endif
elif code == "A"
matrix Rd = b.Rd0
err = add_constraint(&Rd, n, r, c, d)
if !err
b.Rd0 = Rd
endif
elif code == "B"
matrix Rd = b.Rd1
err = add_constraint(&Rd, n, r, c, d)
if !err
b.Rd1 = Rd
endif
elif code == "Adiag"
matrix Rd = b.Rd0
if ok(r)
b.Rd0 = Rd | diag_Rd(n, r)
else
b.Rd0 = Rd | free_diag_Rd(n)
endif
elif code == "Bdiag"
matrix Rd = b.Rd1
if ok(r)
b.Rd1 = Rd | diag_Rd(n, r)
else
b.Rd1 = Rd | free_diag_Rd(n)
endif
endif
## Inform the user if the restriction failed.
/*
At this point it should hold that:
err == 0 : add_constraint worked ok, -- or Adiag/Bdiag: is it conceivable that this happens
together with other A/B restrictions? Then it probably should
also be checked in principle (but doesn't happen here).
*/
if err == -1
printf "Imposing restriction failed, bad input to "
printf "add_constraint.\n"
elif err == 10
printf "Imposing restriction failed, conflicting with "
printf "earlier restrictions.\n"
elif err == 20
printf "Imposing restriction failed, redundant.\n"
endif
if err
printf "(Code %s, ", code
if ok(r)
printf "element %d,%d restricted to %f.)\n", r,c,d
else
printf "no manual restriction.)\n"
endif
endif
return err # 0, -1, 10, or 20
</code>
</gretl-function>
<gretl-function name="SVAR_ident" type="scalar">
<params count="2">
<param name="b" type="bundleref"/>
<param name="verbose" type="int" default="0"/>
</params>
<code>type = b.type
n = b.n
n2 = n*n
ret = 1
if (type==1) || (type==2) || (type == 4) # plain, C-model or SVEC
matrix Rd = get_full_Rd(&b, 0)
matrix Ra = I(n2)
matrix da = vec(I(n))
matrix Rb = Rd[,1:n2]
matrix db = Rd[,n2+1]
elif type == 3 # AB - Model
matrix bRd = b.Rd1 # restrictions on B
matrix aRd = b.Rd0 # restrictions on A
matrix Ra = aRd[,1:n2]
matrix da = aRd[,n2+1]
matrix Rb = bRd[,1:n2]
matrix db = bRd[,n2+1]
endif
if verbose
printf "Constraints in implicit form:\n\n"
loop foreach i Ra da Rb db --quiet
printf "$i:\n%4.0f\n", $i
endloop
endif
ret = ident(Ra, da, Rb, db, verbose)
return ret
</code>
</gretl-function>
<gretl-function name="SVAR_estimate" type="scalar">
<params count="2">
<param name="obj" type="bundleref"/>
<param name="verbosity" type="int" default="1"/>
</params>
<code>/*
this function fills the bundle with the estimated structural
matrices and the covariance matrix of their free elements; it also
calls do_IRF at the end so that the structural VMA is stored into
the bundle
*/
scalar type = obj.type
scalar meth = obj.optmeth
scalar n = obj.n
scalar T = obj.T
matrix vcv
scalar errcode = 0
if type == 4 # do VECM
vecm_est(&obj)
else # estimate ordinary VAR
base_est(&obj)
endif
# grab the instantaneous covariance matrix
matrix Sigma = obj.Sigma
scalar obj.LL0 = VARloglik(obj.T, Sigma)
if type == 1
# plain model: just Cholesky decomposition
matrix C = cholesky(Sigma)
matrix param = vech(C')
# compute the covariance matrix for C
matrix Ss = imp2exp(obj.Rd1)
vcv = coeffVCV(Ss[, 1: cols(Ss)-1], &C)
elif (type==2) || (type == 4)
# C models in a broad sense (including SVEC)
matrix fullRd = get_full_Rd(&obj, verbosity) # this also sets obj.C1 !
# try to set some "sensible" initial values
matrix param = init_C(Sigma, fullRd)
if obj.simann > 0
printf "before simann:\n%12.6f\n", param
matrix dat = vec(Sigma) ~ imp2exp(fullRd)
matrix parcpy = param
zzz = simann(parcpy, "loglik(&parcpy, &dat, -1)", obj.simann)
param = parcpy
printf "after simann:\n%12.6f\n", param
endif
# do estimation; note that vcv is estimated inside "estC"
matrix C = estC(&param, Sigma, fullRd, &vcv, &errcode, meth, verbosity)
elif type == 3
# AB-model
matrix E = obj.E # grab the VAR residuals (needed for initialisation)
matrix bRd = obj.Rd1 # restrictions on B
matrix aRd = obj.Rd0 # restrictions on A
# try to set some "sensible" initial values
# (substitute out call to (former) init_AB)
matrix param = is_standard_AB(aRd, bRd) ? stdAB_init(E, aRd, bRd) : nonstdAB_init(E, aRd, bRd)
# do estimation; note that vcv is estimated inside "estAB"
# (no it actually wasn't, now it is <Sven>)
matrices transfer
matrix C = estAB(&param, Sigma, aRd, bRd, &vcv, &errcode, meth, verbosity, &transfer)
endif # which type
/*
Post-estimation; transfers, long-run matrix, over-id test
*/
if errcode
funcerr "Estimation failed"
else # estimation ran fine
# Copy stuff into the bundle
if type == 3
matrix obj.S1 = transfer[1] # A
matrix obj.S2 = transfer[2] # B
scalar obj.ka = transfer[3]
scalar obj.kb = transfer[4]
endif
/*
We now also store C in the bundle because its computation
was repeated several times elsewhere */
matrix obj.C = C
/* (This was obj.S1 for C models instead of obj.C, don't know why...)
matrix obj.S1 = C */
matrix obj.theta = param
matrix obj.vcv = vcv
# Long-run matrix
# Jan 2018: add the type 4 possibility
if ( (type < 3) && ( rows(obj.Rd1l) || obj.calc_lr ) ) || type == 4
# a plain or C model with long-run constraints, or user switch;
# here we now (Oct 2017) calc and save the long-run matrix
# re-use the C1 matrix from above if possible
matrix C1 = (type == 2 || type == 4) ? obj.C1 : C1mat(obj.VARpar)
matrix obj.lrmat = C1 * obj.C # was: obj.S1
endif
# Indicate that estimation is done (still needed?)
obj.step = 2
# Store IRFs into the bundle
doIRF(&obj)
# store the log-likelihood into the bundle
scalar obj.LL1 = VARloglik(obj.T, obj.Sigma, &C)
# calculate the over-id test in any case (not just for verbosity)
# (C should hopefully be correctly depending on type)
overid = (n * (n+1) / 2 - rows(obj.theta))
if overid > 0
LR = 2 * (obj.LL0 - obj.LL1)
matrix obj.LRoid = {LR; overid; pvalue(X, overid, LR)}
# this is: (stat| dof| pv)
endif
if (verbosity > 0)
SVAR_est_printout(&obj)
endif
endif
return errcode
</code>
</gretl-function>
<gretl-function name="SVAR_cumulate" type="scalar">
<params count="2">
<param name="b" type="bundleref"/>
<param name="nv" type="int"/>
</params>
<code>err = (nv>b.n || nv<1) # was nv<0, but 0 makes no sense (?)
if !err
vn = b.Ynames
printf "Variable %s cumulated\n", vn[nv]
b.cumul = b.cumul | {nv}
b.ncumul = b.ncumul + 1
endif
return err
</code>
</gretl-function>
<gretl-function name="SVAR_boot" type="scalar">
<params count="4">
<param name="obj" type="bundleref"/>
<param name="rep" type="int"/>
<param name="alpha" type="scalar"/>
<param name="quiet" type="bool" default="1"/>
</params>
<code>loop foreach i n k T p type --quiet
scalar $i = obj.$i # copy for convenience, # type = obj.type
endloop
scalar h = obj.horizon
scalar n2 = n*n
# meth = obj.optmeth
matrix m = obj.mu # deterministics
/* --- constraints-related matrices ---------------------*/
matrix Rd1 = obj.Rd1
matrix Rd2 = type==3 ? obj.Rd0 : obj.Rd1l
if type == 4
coint_r = cols(obj.jbeta)
matrix J = zeros(n - coint_r, coint_r) | I(coint_r)
endif
matrix X = obj.X # exogenous variables
matrix start = obj.Y
start = start[1:p,] # Y0
matrix E = cdemean(obj.E) # was: E .- meanc(E) # centre residuals
matrix C = obj.C # (used below at least in maybe_flip_columns(C, &K))
if type == 3 # AB
matrix bmA bmB # needed as memory for transfer
endif
/*
the matrix "bands" will contain the bootstrap results:
each bootstrap replication on one row; each row contains
the vectorisation of the complete IRF matrix
*/
matrix bands = zeros(rep, (h+1) * n2)
matrix U
matrix bmu = (k > 0) ? X*m : {}
bundle bobj = obj # store a copy of the model for bootstrap
obj.nboot = rep # record bootstrap details
obj.boot_alpha = alpha # into original model
if obj.ncumul > 0
matrix to_cum = obj.cumul
matrix tmp = zeros(n, n)
tmp[to_cum,] = 1
sel = selifr(transp(seq(1, n2)), vec(tmp)) # was: n*n
endif
i = 1
failed = 0
set loop_maxiter 16384
# The result matrix (IRFs further below)
if (type > 2) || (rows(obj.Rd1l) || obj.calc_lr )
# save either A,B (for type 3) or C and the long-run matrix
matrix Spar_mat = zeros(rep, 2 * n2)
else # type 1,2,4 and just save the C matrices.
matrix Spar_mat = zeros(rep, n2)
endif
BIASCORR = (type != 4) && obj.biascorr # not available w/unit roots # was: bc
if BIASCORR
matrix Psi = {}
matrix ABC = bias_correction(1000, obj.VARpar, bmu, E, X, start, &Psi) # was: A
endif
printf "\nBootstrapping model (%d iterations)\n\n", rep
flush
loop while i <= rep --quiet
bobj.step = 0 # clear previous bootstrap indicator
/*
generate bootstrap disturbances: first p rows
(corresponding to Y0) are 0; next rows are sampled
with replacement from VAR residuals
*/
U = zeros(p,n) | resample(E)
if rows(bmu) > 0
U = bmu + U
endif
# generate bootstrap data
matrix ABCorA = BIASCORR ? ABC : obj.VARpar # was: A
matrix bY = varsimul(ABCorA, U[p+1:,], start)
matrix bobj.Y = bY # and store it into the bootstrap object
/*
estimate VAR parameters, via VECM if type==4 (KPSW)
or via VAR otherwise
*/
if type == 4
vecm_est(&bobj)
else
base_est(&bobj)
endif
matrix bA = bobj.VARpar # recover estimates (first n rows of companion mat)
matrix bSigma = bobj.Sigma
matrix theta = obj.theta # init original SVAR params (C/A&B apparently, in suitable form...)
# (was: param)
errcode = 0
## Full bias correction
/* (Moved up from below because the new C estimates should be based
on the bc-ed VARpar version.)
But note that the bc-ed VARpar only actually enter the new C matrix if
there are long-run constraints, namely via C1 through fullRd.
Otherwise the new C only depends on Sigma, but we do not update Sigma
in the full biascorr case. (In theory we could, by re-calculating the
residuals using the bc-ed VARpar. We shouldn't, should we?)
*/
if obj.biascorr == 2 && type != 4 # was: bc;
# (The check for SVEC / type!=4 was missing before -- Sven)
scalar H = add_and_smash(&bA, Psi)
# printf("\tH=%g\n",H)
if ok(H)
bobj.VARpar = bA
else
errcode = 101
endif
endif
/* now re-estimate C, according to model type */
if type == 1 # Cholesky
matrix K = cholesky(bSigma)
elif type == 2 || type == 4 # consolidate the SVEC case here
# FIXME: (Sven) The following SVEC check can simply be removed after
# we verify that the call to get_full_Rd() works alright in this case!
# Jan 2018: try to remove it...:
# if type == 4 && rows(bobj.Rd1l)
# funcerr "Bootstrap in SVEC with further long-run restrictions not supported (yet)."
# endif
matrix fullRd = get_full_Rd(&bobj)
matrix K = estC(&theta, bSigma, fullRd, null, &errcode, obj.optmeth, 0)
/* For type==4, this was before:
# FIXME: this is very much WIP
matrix aperp = nullspace(bobj.jalpha')
# set up permanent-transitory constraints
r = cols(bobj.jbeta)
matrix J = zeros(n-r, r) | I(r)
matrix ptRd = (J ** aperp)' ~ zeros(r*(n-r),1)
# matrix J = I(n-r) | zeros(r, n-r)
# moo = (J ** beta)' ~ zeros(r*(n-r),1)
# ptRd |= moo
matrix fullRd = Rd1 | ptRd
*/
elif type == 3 # "AB"
matrices transferAB = array(2) # new: get A,B instead of re-calc'ing it below
matrix K = estAB(&theta, bSigma, Rd2, Rd1, null, &errcode, obj.optmeth, 0, &transferAB)
# elif type == 4 # SVEC, was: "KPSW"
# matrix fullRd = get_full_Rd(&bobj)
# matrix K = estC(&theta, bSigma, fullRd, null, &errcode, obj.optmeth, 0)
endif
## Process and store the simulated C results
if !errcode && rows(K) == n
bobj.step = 2
bobj.theta = theta
# we don't treat the AB-model specially here (there's no reason to)
maybe_flip_columns(C, &K)
if (type == 1) || (type == 2) || (type == 4)
bobj.C = K # was: bobj.S1 = K / is used in doIRF()
Spar_mat[i, 1:n2] = vec(K)' # was just [i,]
/* New Oct 2017: Also bootstrap the long-run matrix if wanted
Jan 2018: add type 4 */
if ( type < 3 && ( rows(bobj.Rd1l) || bobj.calc_lr ) ) || type == 4
# (a plain or C model with long-run constraints, or user switch)
# long-run matrix (C1 comes from get_full_Rd() above (except type 1)):
matrix C1 = (type == 2 || type == 4) ? bobj.C1 : C1mat(bobj.VARpar)
matrix bobj.lrmat = C1 * bobj.C
# attach it to the other bootstrap result
Spar_mat[i, n2+1 : 2*n2] = vec(bobj.lrmat)'
endif
elif type == 3
# (Sven): the following stuff could be gotten from estAB above
bobj.S1 = transferAB[1] # was: bmA
bobj.S2 = transferAB[2] # bmB
Spar_mat[i,] = vec(bobj.S1)' ~ vec(bobj.S2)' # was: vec(bmA)' ~ vec(bmB)'
endif
### This was the place where the full bias correction block lived before ###
endif
if !errcode && rows(K) == n
doIRF(&bobj)
bands[i,] = vec(bobj.IRFs)' # was vec(ir) from ir = bobj.IRFs
i++
else
failed++
outfile stderr --write
printf "Iter %4d failed (error code = %d)\n", i, errcode
outfile --close
endif
endloop
if !quiet
boot_printout(type, n, rep, failed, &Spar_mat)
endif
# quantiles of bootstrapped IRFs used in graphs
matrix locb = quantile(bands, 1 - alpha)
matrix hicb = quantile(bands, alpha)
matrix mdn = quantile(bands, 0.5)
bundle bootdata = null
bootdata.rep = rep # no of replications
bootdata.alpha = alpha # alpha
bootdata.biascorr = obj.biascorr # type of bias correction
matrix bootdata.lo_cb = mshape(locb, h+1, n2) # was: locb # lower bounds
matrix bootdata.hi_cb = mshape(hicb, h+1, n2) # was: hicb # upper bounds
matrix bootdata.mdns = mshape(mdn, h+1, n2) # was: mdn # medians
bundle obj.bootdata = bootdata
return failed
</code>
</gretl-function>
<gretl-function name="SVAR_hd" type="list">
<params count="2">
<param name="Mod" type="bundleref"/>
<param name="nv" type="scalar"/>
</params>
<code>list ret = null
loop foreach i n p t1 t2 type --quiet
scalar $i = Mod.$i
endloop
if nv<0 && nv>n
printf "Hm. There are %d variables in the model. You asked for no. %d\n", n, nv
# (And yno doesn't exist?)
return ret
endif
# compute the exogenous part
if (type < 4)
# this is easy
matrix m = Mod.X * Mod.mu
else
# here we have to take into account the "5 cases"
scalar dcase = Mod.jcase
scalar T = Mod.T
matrix mreg = (dcase == 1) ? {} : ones(T,1)
if dcase > 3
matrix mreg ~= seq(1,T)'
endif
# printf "%8.3f\n", mreg
# printf "%8.3f\n", Mod.X
# printf "%8.3f\n", Mod.mu
matrix m = (mreg ~ Mod.X) * Mod.mu
endif
matrix E = Mod.E
matrix VARpar = Mod.VARpar
if (type == 1) || (type == 2) || (type == 4)
matrix C = Mod.C
elif type == 3
if inbundle(Mod, "C")
matrix C = Mod.C
else
matrix C = Mod.S1 \ Mod.S2
endif
endif
matrix iC = inv(C)
strings Ynames = Mod.Ynames
strings snames = Mod.snames
string yn = Ynames[nv]
smpl t1 t2
if cols(m)>0
Xdet = varsimul(VARpar, m[p+1:,], Mod.Y[1:p,])
else
Xdet = varsimul(VARpar, zeros(Mod.T-p,n), Mod.Y[1:p,])
endif
# printf "nobs = %d, rows(Xdet) = %d\n", $nobs, rows(Xdet)
ret += genseries( sprintf("hd_%s_det", yn), Xdet[,nv])
# the structural shocks
matrix U = E * iC'
rotVARpar = iC * Mod.VARpar * (I(p) ** C)
loop i=1..n --quiet
a = (seq(1,n) .= i)
W = varsimul(rotVARpar, U .* a, zeros(p,n)) * C'
#printf "U[$i]:\n%8.3f\n", W[1:10,]
ret += genseries(sprintf("hd_%s_%s", yn, snames[i]), W[,nv])
endloop
return ret
</code>
</gretl-function>
<gretl-function name="SVAR_coint" type="scalar">
<params count="5">
<param name="SVARobj" type="bundleref"/>
<param name="dcase" type="int" min="1" max="5" default="3"/>
<param name="jbeta" type="matrix"/>
<param name="jalpha" type="matrix"/>
<param name="verbose" type="bool" default="0"/>
</params>
<code>/*
This function doesn't do very much, except setting
up the model for subsequent VECM estimation; "dcase" tells you
which of the "five cases" we want (no constant, restricted
constant, etc), jbeta is simply checked for dimensions and then
copied into the object.
As for jalpha, if it's an empty matrix, that means "just estimate it unrestrictedly", and we set up a
flag accordingly. Otherwise, it's taken to be pre-set to some
fixed value; the intermediate case (contraints on alpha) is not
handled, and I doubt it will ever be.
While we're at it, we also label the structural shocks as
"Perm_1", "Perm_2", "Trans_1", "Trans_2", etc.
*/
scalar n = SVARobj.n
# syntax check
err = (dcase<1) || (dcase>5)
if err
errmsg = sprintf("Invalid dcase value %d", dcase)
funcerr errmsg
endif
# dimensions check
if dcase%2 # nice, huh?
err = err || (rows(jbeta) != n)
else
err = err || (rows(jbeta) != n+1)
endif
if err
errmsg = sprintf("jbeta: has %d rows, should have %d", rows(jbeta), dcase%2 ? n : n+1)
funcerr errmsg
endif
r = cols(jbeta)
err = err || (n < r) # should this be <=? hmm.
# rank check
err = err || (rank(jbeta) < r)
# now check if alpha is ok
d = rows(jalpha)
# d==0 is ok, we'll estimate alpha later
free_a = (d==0)
if !free_a
err = err || (d != n) || (cols(jalpha) != r)
endif
# if anything goes wrong, return
if err
outfile stderr --write --quiet
printf "SVAR_coint: returning on err = %d\n", err
outfile --close
return err
endif
# fill up the object with the info
SVARobj.crank = r
SVARobj.jcase = dcase
SVARobj.jbeta = jbeta
SVARobj.jalpha = jalpha
SVARobj.free_a = free_a
if verbose
if dcase == 1
printf "No constant, "
elif dcase == 2
printf "Restricted constant, "
elif dcase == 3
printf "Unrestricted constant, "
elif dcase == 4
printf "Restricted trend, "
elif dcase == 5
printf "Unrestricted trend, "
endif
printf "beta =\n%9.5f\n", jbeta
if free_a
printf "alpha is unrestricted\n"
else
printf "alpha =\n%9.5f\n", jalpha
printf "PI =\n%9.5f\n", jalpha * jbeta'
endif
endif
# relabel transitory structural shocks
strings sn = SVARobj.snames
if n-r == 1
sn[1] = sprintf("Permanent")
if r == 1
sn[2] = "Transitory"
else
loop i=2..n --quiet
sn[i] = sprintf("Transitory_%d", i-1)
endloop
endif
else
loop i=1..n --quiet
sn[i] = sprintf("Permanent_%d", i)
endloop
if r == 1
sn[n] = "Transitory"
else
loop i=1..r --quiet
sn[n-r+i] = sprintf("Transitory_%d", i)
endloop
endif
endif
SVARobj.snames = sn
return err
</code>
</gretl-function>
<gretl-function name="GetShock" type="series">
<params count="2">
<param name="SVARobj" type="bundleref"/>
<param name="i" type="scalar"/>
</params>
<code>series ret = NA
scalar n = SVARobj.n
if (i <= n) && (i > 0)
scalar type = SVARobj.type
if (type == 1) || (type == 2) || (type == 4)
matrix C = SVARobj.C # was: SVARobj.S1
elif type == 3
if inbundle(SVARobj, "C") # maybe not yet computed
matrix C = SVARobj.C
else
matrix C = SVARobj.S1 \ SVARobj.S2
endif
endif
matrix iC = inv(C')
scalar extra = $nobs - rows(SVARobj.E)
if extra > 0
set warnings off
matrix tmp = ones(extra,1) .* NA
else
matrix tmp = {}
endif
tmp |= SVARobj.E * iC[,i]
series ret = tmp
snames = SVARobj.snames # strings array?
string vlab = snames[i]
setinfo ret --description="@vlab"
endif
return ret
</code>
</gretl-function>
<gretl-function name="IRFplot" type="void">
<params count="4">
<param name="obj" type="bundleref"/>
<param name="snum" type="int"/>
<param name="vnum" type="int"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
</params>
<code>IRFsave("display", &obj, snum, vnum, keypos)
</code>
</gretl-function>
<gretl-function name="IRFsave" type="void">
<params count="5">
<param name="outfilename" type="string"/>
<param name="obj" type="bundleref"/>
<param name="snum" type="int"/>
<param name="vnum" type="int"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
</params>
<code># negative snum is allowed and means to flip the shock
# copy and prepare some input
string tmpfile, tmpout
n = obj.n
matrix IRFmat = obj.IRFs
scale = 1 # possibly changed later
if obj.normalize == 1
matrix tmp = mshape(IRFmat[1,], n, n)
endif
if obj.nboot
boot = obj.bootdata # bundle?
bc = boot.biascorr
endif
scalar sfrom sto vfrom vto
is_srange = range(snum, n, &sfrom, &sto)
is_vrange = range(vnum, n, &vfrom, &vto)
# cycle through all (selected) shocks and variables
loop snum = sfrom..sto -q
# do checks
err = check_bounds(snum, 1, n)
if err == 1
printf "Shock number %d out of bounds\n", abs(snum) # abs because of flipping
return
endif
string sn = obj.snames[abs(snum)]
# normalization / scaling
if obj.normalize == 1
scale = tmp[abs(snum), abs(snum)]
endif
loop vnum = vfrom..vto -q
# do checks
err = check_bounds(1, vnum, n)
if err == 2
printf "Variable number %d out of bounds\n", vnum
return
endif
string yn = obj.Ynames[vnum]
# produce plots
if obj.ncumul == 0
if obj.nboot == 0
tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos)
else
tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, null, &boot, bc)
endif
else
matrix cumul = obj.cumul
if obj.nboot == 0
tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, &cumul)
else
tmpfile = IRFgrph(IRFmat, snum, vnum, scale, sn, yn, keypos, &cumul, &boot, bc)
endif
endif
if (outfilename == "display") || (sfrom == sto && vfrom == vto) # single plot, no indices
tmpout = outfilename
else
strings be = basename(outfilename) # gives 2-elem array
if is_vrange # several v indices
if is_srange # several s indices
tmpout = sprintf("%s_%d%d.%s", be[1], snum, vnum, be[2])
else
tmpout = sprintf("%s_%d.%s", be[1], vnum, be[2])
endif
elif is_srange # several s indices
tmpout = sprintf("%s_%d.%s", be[1], snum, be[2])
else
funcerr "shouldn't happen"
endif
endif
gnuplot --input="@tmpfile" --output="@tmpout"
endloop
endloop
</code>
</gretl-function>
<gretl-function name="FEVD" type="matrix">
<params count="1">
<param name="SVARobj" type="bundleref"/>
</params>
<code>n = SVARobj.n
h = SVARobj.horizon + 1
sIRF = SVARobj.IRFs
matrix ret = zeros(h, n*n)
ctmp = cum(sIRF .* sIRF)
loop i=1..h --quiet
tmp = mshape(ctmp[i,],n,n)'
ret[i,] = vec(tmp ./ sumc(tmp))'
endloop
return ret
</code>
</gretl-function>
<gretl-function name="GUI_SVAR" type="bundle" pkg-role="gui-main">
<params count="16">
<param name="type" type="int" min="1" max="3" default="1">
<description>Model type</description>
<labels count="3">
"plain (Cholesky)" "C-model" "AB-model" </labels>
</param>
<param name="Y" type="list">
<description>VAR variables</description>
</param>
<param name="X" type="list" optional="true">
<description>Exogenous regressors</description>
</param>
<param name="hasconst" type="bool" default="1">
<description>Constant</description>
</param>
<param name="hastrend" type="bool" default="0">
<description>Time trend</description>
</param>
<param name="hasseas" type="bool" default="0">
<description>Seasonal dummies</description>
</param>
<param name="l" type="int" min="1" default="1">
<description>Lags</description>
</param>
<param name="h" type="int" min="0">
<description>Horizon</description>
</param>
<param name="R1" type="matrixref" optional="true">
<description>Restriction pattern (short-run C or B)</description>
</param>
<param name="R2" type="matrixref" optional="true">
<description>Restriction pattern (long-run C or A)</description>
</param>
<param name="b" type="int" min="0">
<description>Bootstrap replications</description>
</param>
<param name="alpha" type="scalar" min="0" max="1" default="0.9">
<description>Bootstrap alpha</description>
</param>
<param name="biascorr" type="int" min="0" max="2" default="0">
<description>Bias correction</description>
<labels count="3">
"None" "Partial" "Full" </labels>
</param>
<param name="checkident" type="bool" default="0">
<description>Check identification</description>
</param>
<param name="cumix" type="matrixref" optional="true">
<description>Indices of responses to cumulate</description>
</param>
<param name="optmeth" type="int" min="0" max="4" default="4">
<description>Optimization method</description>
<labels count="5">
"BFGS (numerical score)" "BFGS (analytical score)" "Newton-Raphson (numerical score)" "Newton-Raphson (analytical score)" "Scoring algorithm" </labels>
</param>
</params>
<code>n = nelem(Y)
# stick together deterministics and other exog.
list lX = dropcoll(determ(X, hasconst, hastrend, hasseas))
# initialize the model bundle
bundle m = SVAR_setup(modelstring(type), Y, lX, l)
if h > 0
m.horizon = h
endif
# copy options (new by Sven)
m.biascorr = biascorr
m.optmeth = optmeth
## implement the cumulation spec (sven 1.0.2)
if !isnull(cumix)
# ensure column vector
cumix = vec(cumix)
# input checks (numbers out of bounds)
if max(cumix) > nelem(Y) || min(cumix) < 1
print "Invalid cumulation specification!"
print "(No responses will be cumulated.)"
else # sensible cumulation spec
loop i=1..rows(cumix) -q
SVAR_cumulate(&m, cumix[i])
endloop
endif
endif
## process restrictions
if type == 1
if !isnull(R1) || !isnull(R2)
print "Estimating plain model. Discarding provided restrictions."
endif
if checkident
print "(Identification trivially given in plain model.)"
endif
else # C or AB model
# input check
if isnull(R1) && isnull(R2)
funcerr "Must provide some restrictions for C and AB models!"
endif
if type == 3 && ( isnull(R1) || isnull(R2) )
funcerr "Must provide restrictions on A and B for AB model!"
endif
# transform the R1-matrix to SVAR-style restrictions
if !isnull(R1)
r = rows(R1)
c = cols(R1)
if (r != n || c != n)
funcerr "wrong R1 dimensions"
endif
string sBorC = type == 3 ? "B" : "C" # new by Sven 1.0.2
loop i=1..n -q
loop j=1..n -q
scalar rij = R1[i,j]
if ok(rij) # valid number = restricted element
SVAR_restrict(&m, sBorC, i, j, rij)
endif
endloop
endloop
endif
# still need to consider the A or longrun-C matrix
# transform the R2-matrix to SVAR-style restrictions
if !isnull(R2)
r2 = rows(R2)
c2 = cols(R2)
if (r2 != n || c2 != n)
funcerr "wrong R2 dimension"
endif
string sAorlrC = type == 3 ? "A" : "lrC"
loop i=1..n -q
loop j=1..n -q
scalar rij = R2[i,j]
if ok(rij) # valid number = restricted element
SVAR_restrict(&m, sAorlrC, i, j, rij)
endif
endloop
endloop
endif
endif
# do an explicit ID check (sven 1.0.2)
scalar id_ok = 1
if checkident
if (type == 2 && !isnull(R2) ) # longrun C
print "FIXME: not yet implemented for models with long-run restrictions"
else
print "Check identification:"
id_ok = SVAR_ident(&m, 1) # request verbosity==1 to get messages
endif
endif
# and of course estimate
if !id_ok
return m
else
SVAR_estimate(&m)
if b > 0
SVAR_boot(&m, b, alpha, 0)
endif
endif
return m
</code>
</gretl-function>
<gretl-function name="GUI_plot" type="void" pkg-role="bundle-plot">
<params count="2">
<param name="b" type="bundleref"/>
<param name="ptype" type="int" min="0" max="2" default="0">
<description>Plot type</description>
<labels count="3">
"IRF" "FEVD" "Historical decomposition" </labels>
</param>
</params>
<code>scalar n = b.n
matrix tt = seq(0, b.horizon)'
if ptype == 0
IRFplot(&b, 0, 0) # all in one go
elif ptype == 1
FEVDplot(&b, 0)
elif ptype == 2
HDplot(&b, 0)
endif
</code>
</gretl-function>
<gretl-function name="FEVDplot" type="void">
<params count="3">
<param name="obj" type="bundleref"/>
<param name="vnum" type="int" default="0"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
</params>
<code>FEVDsave("display", &obj, vnum, keypos)
</code>
</gretl-function>
<gretl-function name="FEVDsave" type="void">
<params count="4">
<param name="outfilename" type="string"/>
<param name="obj" type="bundleref"/>
<param name="vnum" type="int" default="0"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
</params>
<code>scalar n = obj.n
scalar vfrom vto
is_vrange = range(vnum, n, &vfrom, &vto)
matrix Fmat = FEVD(&obj)
loop vnum = vfrom..vto -q
err = check_bounds(1, vnum, n)
if err == 2
printf "Variable number %d out of bounds\n", vnum
return
endif
string tmpout = (outfilename == "display") ? "display" : sprintf("%s_%d", outfilename, vnum)
string tmpfile = FEVDgrph(Fmat, vnum, obj.Ynames[vnum], obj.snames, keypos)
gnuplot --input="@tmpfile" --output="@tmpout"
endloop
</code>
</gretl-function>
<gretl-function name="HDplot" type="void">
<params count="2">
<param name="obj" type="bundleref"/>
<param name="vnum" type="int" default="0"/>
</params>
<code>HDsave("display", &obj, vnum)
</code>
</gretl-function>
<gretl-function name="HDsave" type="void">
<params count="3">
<param name="outfilename" type="string"/>
<param name="obj" type="bundleref"/>
<param name="vnum" type="int" default="0"/>
</params>
<code># more fashionable (=Dynare-like) style
# Interpret vnum==0 as meaning "all 1..n"
scalar n = obj.n
scalar vfrom vto
is_vrange = range(vnum, n, &vfrom, &vto)
loop vnum = vfrom..vto -q
err = check_bounds(1, vnum, n)
if err == 2
printf "Variable number %d out of bounds\n", vnum
return
endif
list tmp = SVAR_hd(&obj, vnum)
tmp -= tmp[1] # take away deterministic component
loop i=1..nelem(tmp) --quiet
string sn = obj.snames[i]
setinfo tmp[i] --graph-name="@sn"
endloop
series tmpvar = sum(tmp)
string gnam = sprintf("%s (stoch. component)", obj.Ynames[vnum])
setinfo tmpvar --graph-name="@gnam"
tmp = tmpvar tmp # put the line with the target var name first
string tmpout = (outfilename == "display") ? "display" : sprintf("%s_%d", outfilename, vnum)
plot tmp
option time-series
option single-yaxis
option with-boxes
option with-lines=tmpvar
printf "set title \"HD for %s\"", obj.Ynames[vnum]
end plot --output="@tmpout"
endloop
</code>
</gretl-function>
<gretl-function name="SVAR_bundle_print" type="void" pkg-role="bundle-print">
<params count="1">
<param name="b" type="bundleref"/>
</params>
<code># Some specification echoing (sven 1.0.2)
# (not sure whether this should go into SVAR_est_printout() instead...)
loop foreach i type n k p --quiet
scalar $i = b.$i
endloop
printf "Model type: %s\n", modelstring(type)
strings Ynames = b.Ynames
print "Endogenous variables:"
loop i=1..n --quiet
printf "%s", Ynames[i]
if (i==n)
printf "\n"
else
printf ", "
endif
endloop
printf "\n"
if k>0
print "Exogenous variables:"
strings Xnames = b.Xnames
loop i=1..k --quiet
printf "%s", Xnames[i]
if (i==k)
printf "\n"
else
printf ", "
endif
endloop
else
print "(No exogenous variables.)"
endif
printf "\n"
printf "Restriction patterns:\n\n"
if type == 1
print "Lower-triangular C matrix (Choleski decomposition)\n"
elif type == 2 # C-model
if rows(b.Rd1)>0
printf "Short-run restrictions:\n"
printf "%3.0f\n", b.Rd1
else
print "(No short-run restrictions)"
endif
if rows(b.Rd1l)>0
print "Long-run restrictions:"
printf "%3.0f\n", b.Rd1l
else
print "(No long-run restrictions.)"
endif
elif type == 3 # AB-model
if rows(b.Rd0)
print "Restrictions on A:"
printf "%3.0f\n", b.Rd0
else
print "(No restrictions on A.)"
endif
if rows(b.Rd1)
print "Restrictions on B:"
printf "%3.0f\n", b.Rd1
else
print "(No restrictions on B.)"
endif
endif
# only print this if actual estimation took place
if inbundle(b, "Sigma")
printf "Sigma = \n%10.6f", b.Sigma
SVAR_est_printout(&b)
endif
</code>
</gretl-function>
<gretl-function name="determ" type="list" private="1">
<params count="4">
<param name="X" type="list"/>
<param name="cnst" type="bool"/>
<param name="trnd" type="bool"/>
<param name="seas" type="bool"/>
</params>
<code>list ret = cnst ? const : null
if trnd
list ret += time
endif
if seas
ipd = 1/$pd
tt = time % $pd
loop i=1..$pd-1 --quiet
ret += genseries(sprintf("seas_%d", i), (tt == i) - ipd)
endloop
endif
# stick together deterministics and other exog.
list ret = ret || X
return ret
</code>
</gretl-function>
<gretl-function name="basename" type="strings" private="1">
<params count="1">
<param name="fn" type="string"/>
</params>
<code>string base = regsub(fn, "(.*)\.([^\.]*)", "\1")
string ext = fn + (strlen(base) + 1)
return defarray(base, ext)
</code>
</gretl-function>
<gretl-function name="range" type="scalar" private="1">
<params count="4">
<param name="a" type="scalar"/>
<param name="n" type="scalar"/>
<param name="n0" type="scalarref"/>
<param name="n1" type="scalarref"/>
</params>
<code># this is used in the plotting function to have
# 0 as a synonym for "all"
#
# The flipping of the shock by passing a negative number
# is then not possible to specify; users have to do it explicitly
# in their own loop then.
if a != 0
ret = 0
n0 = a
n1 = a
else
ret = 1
n0 = 1
n1 = n
endif
return ret
</code>
</gretl-function>
<gretl-function name="base_est" type="scalar" private="1">
<params count="1">
<param name="SVARobj" type="bundleref"/>
</params>
<code>/*
takes a SVAR object and adds the basic VAR estimates
to the bundle; returns an errcode (virginity check).
*/
scalar step = SVARobj.step
err = (step>0)
if err
printf "Base estimation done already!\n"
else
matrix mY = SVARobj.Y
scalar p = SVARobj.p
scalar n = SVARobj.n
scalar k = SVARobj.k
scalar T = SVARobj.T
matrix E
matrix mreg = SVARobj.X ~ mlag(mY, seq(1,p))
matrix olspar = mols(mY[p+1:T,], mreg[p+1:T,], &E)
matrix Sigma = (E'E) ./ (T - (n*p + k)) # was: / df
matrix SVARobj.VARpar = transp(olspar[k+1:k+n*p,])
matrix SVARobj.mu = (k>0) ? olspar[1:k,] : {} # was: d
matrix SVARobj.Sigma = Sigma
matrix SVARobj.E = E
SVARobj.step = 1
SVARobj.LL0 = T * (n* log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
endif
return err
</code>
</gretl-function>
<gretl-function name="vecm_det" type="matrix" private="1">
<params count="2">
<param name="T" type="int"/>
<param name="dcase" type="int"/>
</params>
<code># build the deterministic matrix for the VECM; if alpha is
# empty, it will be estimated via ols later
# deterministics
# note that in the "even cases" (restr. const or trend)
# the restricted term is _not_ included in x, since its
# parameter will be recovered later via alpha*beta
matrix mreg = (dcase < 3) ? {} : ones(T,1)
if dcase == 5
matrix mreg ~= seq(1,T)'
endif
return mreg
</code>
</gretl-function>
<gretl-function name="vecm_est" type="scalar" private="1">
<params count="1">
<param name="SVARobj" type="bundleref"/>
</params>
<code>/*
We can't afford to be too flexible here, and the intended usage
is as follows.
We assume the model has already been declared as a SVEC (type=4)
model. We also assume that the cointegration properties (beta mat plus
deterministics) have already been set up via the SVAR_coint() function, so that we already have beta and alpha available as "jbeta" and
"jalpha", respectively. Finally, we take care of proper treatment of
deterministics, via the scalar "jcase" (1 to 5).
*/
# --- preliminary checks
err = (SVARobj.step > 0)
if err
printf "Base estimation done already!\n"
return err
endif
err = (SVARobj.type != 4)
if err
printf "Wrong model type!\n"
return err
endif
# --- grab stuff from the bundle
matrix mY = SVARobj.Y
matrix jbeta = SVARobj.jbeta
matrix jalpha = SVARobj.jalpha
scalar p = SVARobj.p
scalar n = SVARobj.n
scalar k = SVARobj.k
scalar r = cols(jbeta)
scalar dcase = SVARobj.jcase
scalar ols_al = rows(jalpha) == 0
scalar T = SVARobj.T
# --- first thing: do we have a pre-set value for alpha?
matrix dY = diff(mY)
matrix dep = dY[p+1:T,]
matrix E = {}
# deterministics
matrix mreg = vecm_det(T, dcase)
# ECM terms
matrix ECM = mlag(mY * jbeta[1:n,],1)
if dcase == 2
ECM = ECM .+ jbeta[n+1, ]
elif dcase == 4
ECM += seq(1,T)'jbeta[n+1, ]
endif
if ols_al
# alpha must be estimated together with the rest of the VECM
matrix mreg ~= SVARobj.X ~ ECM
else
matrix dep -= (ECM[p+1:T,] * jalpha')
matrix mreg ~= SVARobj.X
endif
# extra lags
if p > 1
mreg ~= mlag(dY, seq(1,p-1))
endif
# trim the first rows
if rows(mreg)>0
mreg = mreg[p+1:T,]
# printf "%d rows\n", rows(dep)
# printf "%8.3f\n", (dep[1:10,] ~ mreg[1:10,])
matrix olspar = mols(dep, mreg, &E)
else
matrix olspar = {}
E = dep
endif
df = T - (n*p + k - (n-r))
matrix Sigma = (E'E)./T
# print olspar Sigma
# --- construct the various matrices required later
rp = rows(olspar)
# alpha first (should be before the gammas, if estimated)
if ols_al
ng = n*(p-1)
jalpha = olspar[rp-ng-r+1:rp-ng,]'
endif
# exogenous
if dcase == 1
matrix mu = {}
scalar nd = 0
elif dcase == 2
matrix mu = jbeta[n+1,] * jalpha' # was: jalpha * jbeta[n+1,] # FIXME
scalar nd = 0
elif dcase == 3
matrix mu = olspar[1,]
scalar nd = 1
elif dcase == 4
matrix mu = olspar[1,] | (jbeta[n+1,] * jalpha') # was: (jalpha * jbeta[n+1,])' # FIXME
scalar nd = 1
elif dcase == 5
matrix mu = olspar[1:2,]
scalar nd = 2
endif
if k > 0
matrix sel = nd + seq(0, k-1)
mu = mu ~ olspar[sel,]
endif
/*
companion form in levels
*/
Pi = jalpha * jbeta[1:n,]'
if p>1
if ols_al
# the Gammas should be after the ECM terms
ini = nd + r + 1
else
# the Gammas should be right after the deterministics
ini = nd + 1
endif
fin = ini + (p-1) * n - 1
matrix A = olspar[ini:fin,] | zeros(n,n)
A += I(n) + Pi' | -olspar[ini:fin,]
else
matrix A = I(n) + Pi'
endif
matrix SVARobj.VARpar = A'
matrix SVARobj.mu = mu
matrix SVARobj.Sigma = Sigma
matrix SVARobj.E = E
matrix SVARobj.jalpha = jalpha
SVARobj.step = 1
SVARobj.LL0 = T*(n*log(2*$pi) - 0.5*log(det(Sigma)) - 0.5)
return err
</code>
</gretl-function>
<gretl-function name="N2_ify" type="matrix" private="1">
<params count="1">
<param name="A" type="matrix"/>
</params>
<code>n2 = rows(A)
n = int(sqrt(n2))
k = cols(A)
matrix e = vec(transp(mshape(seq(1,n2),n,n)))
return A + A[e,1:k]
</code>
</gretl-function>
<gretl-function name="has_unit_diag" type="scalar" private="1">
<params count="1">
<param name="Rd" type="matrix"/>
</params>
<code>ret = 0
n2 = cols(Rd) - 1
n = sqrt(n2)
matrix test = transp(seq(0,n2-1)) .= seq(0,n2-1,n+1)
test |= ones(1,n)
matrix e
mols(test, Rd', &e)
# ret = maxr(sumc(e.*e)) < 1.0e-12
return maxr(sumc(e.*e)) < 1.0e-12 # ret
</code>
</gretl-function>
<gretl-function name="has_free_diag" type="scalar" private="1">
<params count="1">
<param name="Rd" type="matrix"/>
</params>
<code>n2 = cols(Rd) - 1
n = sqrt(n2)
matrix e = 1 + seq(0,n2-1,n+1)
matrix test = Rd[,e]
# ret = sumc(abs(test)) < 1.0e-12
return ( sumc(abs(test)) < 1.0e-12 ) # ret
</code>
</gretl-function>
<gretl-function name="mat_exp" type="matrix" private="1">
<params count="3">
<param name="theta" type="matrix"/>
<param name="Ss" type="matrix"/>
<param name="force_pos" type="bool" default="0"/>
</params>
<code># FIXME: force_pos switch is not used currently (commented out below)
# we don't check for conformability, but
# cols(Ss) should be equal to rows(theta)+1
n2 = rows(Ss)
n = int(sqrt(n2))
k = cols(Ss) - 1
matrix C = (k>0) ? ( Ss[,1:k]*theta + Ss[,k+1] ) : Ss
C = mshape(C,n,n)
/*
if force_pos
neg = (diag(C)') .< 0
pos = !neg
C = C .* (pos - neg)
endif
*/
return C
</code>
</gretl-function>
<gretl-function name="unscramble_dat" type="void" private="1">
<params count="3">
<param name="dat" type="matrixref"/>
<param name="Sigma" type="matrixref"/>
<param name="Ss" type="matrixref"/>
</params>
<code>n2 = rows(dat)
n = int(sqrt(n2))
k = cols(dat)
Sigma = mshape(dat[,1],n,n)
Ss = dat[,2:k]
</code>
</gretl-function>
<gretl-function name="VARloglik" type="scalar" private="1">
<params count="3">
<param name="T" type="scalar"/>
<param name="Sigma" type="matrix"/>
<param name="C" type="matrixref" optional="true"/>
</params>
<code># computes the (concentrated) loglikelihood for a VAR
# the matrix C (such that in theory CC' = Sigma) may be null,
# for just-identified models
n = rows(Sigma)
ll = n * 1.83787706640934548355 # ln(2*$pi)
if isnull(C) # just-identified model
ll += log(det(Sigma)) + n
else
matrix KK = invpd(C*C')
ll += -log(det(KK)) + tr(KK*Sigma)
endif
return -(T/2) * ll
</code>
</gretl-function>
<gretl-function name="loglik" type="scalar" private="1">
<params count="3">
<param name="theta" type="matrixref"/>
<param name="dat" type="matrixref"/>
<param name="modeltype" type="scalar"/>
</params>
<code># modeltype < 0->C; modeltype>=0: AB (contains free elements of a)
matrix Sigma = {}
matrix Ss = {}
unscramble_dat(&dat, &Sigma, &Ss)
if modeltype < 0
matrix C = mat_exp(theta, Ss, 0)
else
p1 = modeltype
p2 = rows(theta)
matrix aSs = Ss[,1:p1+1]
matrix bSs = Ss[,p1+2:p2+2]
matrix A B
ABmat_exp(theta, aSs, bSs, &A, &B)
# was: matrix C = B/A
matrix C = A\B
endif
matrix KK = invpd(C*C')
ll = det(KK) # should always be positive
ll = (ll<=0) ? NA : -0.5 * (tr(KK*Sigma) - log(ll))
return ll
</code>
</gretl-function>
<gretl-function name="InfoMat" type="matrix" private="1">
<params count="3">
<param name="CorB" type="matrix"/>
<param name="S" type="matrix"/>
<param name="A" type="matrixref" optional="true"/>
</params>
<code>/*
merged from InfoMatC and InfoMatAB
First case C model (A is null), latter AB model.
*/
matrix tmp = isnull(A) ? I(rows(CorB)) : (A\CorB) | -I(rows(A))
tmp = S' (tmp ** inv(CorB'))
return tmp * N2_ify(tmp')
</code>
</gretl-function>
<gretl-function name="coeffVCV" type="matrix" private="1">
<params count="3">
<param name="S" type="matrix"/>
<param name="rhs" type="matrixref"/>
<param name="lhs" type="matrixref" optional="true"/>
</params>
<code># C or AB model
matrix IM = isnull(lhs) ? InfoMat(rhs, S) : InfoMat(rhs, S, &lhs)
# (should be correct with new InfoMat)
# quick-dirty check for singularity
if rcond(IM)>1.0e-10
matrix iIM = invpd(IM)/$nobs
else
matrix evec
l = eigensym(IM, &evec)
printf "\n\nInformation matrix is not pd!!\n\n%12.5f\n", IM
printf "Eigenvalues:\n%12.5f\n", l
matrix e = (l .< 1.0e-07)
printf "Troublesome eigenvectors:\n%12.5f\n", selifc(evec, e')
printf "S:\n%4.0f\n", S
matrix iIM = zeros(rows(IM), cols(IM))
endif
return qform(S,iIM)
</code>
</gretl-function>
<gretl-function name="maybe_flip_columns" type="void" private="1">
<params count="2">
<param name="C" type="matrix"/>
<param name="X" type="matrixref"/>
</params>
<code>/*
the objective here is to make X as similar as possible to C
by flipping the sign of the columns. Used for bootstrapping, to have IRFs with comparable signs.
*/
n = rows(C)
matrix sel = seq(1,n)
matrix plus = sumc((C + X).^2)
matrix minus = sumc((C - X).^2)
matrix flip = plus .< minus
if sumr(flip) > 0
sel = selifc(sel, flip)
X[,sel] = -X[,sel]
endif
</code>
</gretl-function>
<gretl-function name="printStrMat" type="void" private="1">
<params count="3">
<param name="X" type="matrix"/>
<param name="V" type="matrix"/>
<param name="name" type="string"/>
</params>
<code>n = rows(X)
matrix x = vec(X)
matrix cfse = vec(X)
matrix se = sqrt(diag(V))
matrix numzero = selifr(seq(1,rows(se))', (se .< 1.0e-15))
if rows(numzero) > 1
se[numzero] = 0.0
endif
cfse ~= se
string parnames = ""
loop j=1..n --quiet
loop i=1..n --quiet
sprintf tmpstr "%s[%2d;%2d]", name, i, j
parnames += tmpstr
if (j<n) || (i<n)
parnames += ","
endif
endloop
endloop
modprint cfse parnames
</code>
</gretl-function>
<gretl-function name="SVAR_est_printout" type="void" private="1">
<params count="1">
<param name="mod" type="bundleref"/>
</params>
<code># scalar type = mod.type
printf "Optimization method = "
strings os = defarray("BFGS (numerical)", "BFGS (analytical)", "Newton-Raphson (numerical)", "Newton-Raphson (analytical score)", "Scoring algorithm")
printf "%s\n", os[mod.optmeth + 1] # Because optmeth starts at 0 here
printf "Unconstrained Sigma:\n%12.5f\n", mod.Sigma
if mod.type < 3 || mod.type == 4 # plain or C-model, or (Jan 2018) SVEC
# Standard C matrix
printStrMat(mod.C, mod.vcv, "C") # was mod.S1
# Long-run matrix
if rows(mod.Rd1l) || mod.calc_lr || mod.type == 4
# either long-run restr., or user wish, or SVEC
# (this could be much embellished, probably)
printf "Estimated long-run matrix"
if rows(mod.Rd1l)
printf " (restricted)"
endif
printf "\n"
matrix longrun = mod.lrmat
# round to zero for printout (also do it in general?)
longrun = (abs(longrun) .< 1e-15) ? 0.0 : longrun
print longrun
endif
elif mod.type == 3 # AB, new <Sven>
n2 = (mod.n)^2
if mod.ka > 0
printStrMat(mod.S1, mod.vcv[1:n2, 1:n2], "A")
endif
if mod.kb > 0
printStrMat(mod.S2, mod.vcv[n2+1 : 2*n2, n2+1 : 2*n2], "B")
endif
endif
printf " Log-likelihood = %g\n", mod.LL1
if inbundle(mod, "LRoid") # over-id test was done
printf " Overidentification LR test = %g (%d df, pval = %g)\n\n", mod.LRoid[1], mod.LRoid[2], mod.LRoid[3]
endif
</code>
</gretl-function>
<gretl-function name="add_and_smash" type="scalar" private="1">
<params count="2">
<param name="A" type="matrixref"/>
<param name="Psi" type="matrix" const="true"/>
</params>
<code>matrix Ab = A + Psi
# now check stationarity
scalar maxmod = max_eval(Ab)
h = 0.99
H = 1
maxiter = 1000
iter = 0
loop while (maxmod > 0.9999) && (iter < maxiter) --quiet
iter++
Ab = A + H .* Psi
maxmod = max_eval(Ab)
H *= h
# printf "H = %g\n", H
endloop
A = Ab
H = (iter >= maxiter) ? NA : H
return H
</code>
</gretl-function>
<gretl-function name="is_standard_AB" type="scalar" private="1">
<params count="2">
<param name="aRd" type="matrix"/>
<param name="bRd" type="matrix"/>
</params>
<code>test = has_unit_diag(aRd) && has_free_diag(bRd)
return test
</code>
</gretl-function>
<gretl-function name="scoreAB" type="void" private="1">
<params count="4">
<param name="ret" type="matrixref"/>
<param name="theta" type="matrixref"/>
<param name="dat" type="matrixref"/>
<param name="p1" type="int"/>
</params>
<code># FIXME: This func might benefit from getting the pre-computed C matrix
matrix Sigma Ss A B
unscramble_dat(&dat, &Sigma, &Ss)
p2 = rows(theta)
matrix aSs = Ss[,1:p1+1]
matrix bSs = Ss[,p1+2:p2+2]
n = rows(Sigma)
ABmat_exp(theta, aSs, bSs, &A, &B)
matrix iA = inv(A)
matrix C = iA*B
if p1>0
matrix S = aSs[,1:p1]
S = -((C') ** iA) * S
else
matrix S = {}
endif
if p2>p1
S ~= bSs[,1:cols(bSs)-1]
endif
matrix iC = inv(C)
matrix ret = iC' (qform(iC,Sigma) - I(n))
ret = vec(ret)'S
</code>
</gretl-function>
<gretl-function name="stdAB_init" type="matrix" private="1">
<params count="3">
<param name="mX" type="matrix"/>
<param name="aRd" type="matrix"/>
<param name="bRd" type="matrix"/>
</params>
<code># mX should contain the VAR residuals
matrix aSs = imp2exp(aRd)
matrix bSs = imp2exp(bRd)
n = cols(mX)
p = cols(aSs) - 1
if p > 0
matrix S = aSs[,1:p]
matrix s = aSs[,p+1]
matrix sel = vec(transp(mshape(seq(1,n*n),n,n)))
matrix dep = vec(mX)
matrix reg = (I(n) ** mX)*S[sel,]
matrix e
matrix b = -mols(dep, reg, &e)
e = mshape(e, rows(mX), n)
e = sqrt(meanc(e.^2))'
S = bSs[,1:cols(bSs)-1]
sel = 1 + seq(0,n-1)*(n+1)
matrix a = zeros(n*n,1)
a[sel] = e
e = S'a
matrix ret = b | e
else
matrix Sigma = (mX'mX) / rows(mX)
matrix ret = init_C(Sigma, bRd)
endif
return ret
</code>
</gretl-function>
<gretl-function name="nonstdAB_init" type="matrix" private="1">
<params count="3">
<param name="E" type="matrix"/>
<param name="aRd" type="matrix"/>
<param name="bRd" type="matrix"/>
</params>
<code>matrix aSs = imp2exp(aRd)
matrix bSs = imp2exp(bRd)
T = rows(E)
n = cols(E)
matrix Sigma = (E'E)./T
matrix C = cholesky(Sigma)
startA = (rows(aRd) > rows(bRd))
#startA = 1 # provisional
ka = cols(aSs) - 1
if ka>0
matrix Sa = aSs[,1:ka]
endif
matrix sa = aSs[,ka+1]
kb = cols(bSs) - 1
if kb>0
matrix Sb = bSs[,1:kb]
endif
matrix sb = bSs[,kb+1]
if startA == 1
matrix giSb = invpd(Sb'Sb)*Sb'
matrix tmp = giSb*(C' ** I(n))
if ka>0
matrix giSa = invpd(Sa'Sa)*Sa'
matrix gama = 0.1*ones(ka,1)
matrix gamb = tmp * (Sa*gama + sa) - giSb*sb
tmp = giSa*(inv(C)' ** I(n))
gama = tmp * (Sb*gamb + sb) - giSa*sa
else
matrix gama = {}
matrix gamb = tmp*sa - Sb'sb
endif
else
matrix giSa = invpd(Sa'Sa)*Sa'
matrix tmp = giSa*(inv(C)' ** I(n))
if kb>0
matrix giSb = invpd(Sb'Sb)*Sb'
matrix gamb = 0.1*ones(kb,1)
matrix gama = tmp * (Sb*gamb + sb) - giSa*sa
tmp = giSb*(C' ** I(n))
gamb = tmp * (Sa*gama + sa) - giSb*sb
else
matrix gama = tmp*sb - Sa'sa
matrix gamb = {}
endif
endif
if kb > 0
scale = tr(Sigma) / tr(C*C')
printf "nonstdAB_init: Scale = %g\n", scale
gamb = sqrt(scale) .* gamb
/*
elif ka>0
scalar scale = tr(Sigma) / tr(C*C')
printf "nonstdAB_init: Scale = %g\n", scale
gama = sqrt(scale) ./ gama
*/
endif
matrix ret = gama | gamb
return ret
</code>
</gretl-function>
<gretl-function name="ABmat_exp" type="void" private="1">
<params count="5">
<param name="theta" type="matrix"/>
<param name="aSs" type="matrix"/>
<param name="bSs" type="matrix"/>
<param name="A" type="matrixref"/>
<param name="B" type="matrixref"/>
</params>
<code>scalar p1 = cols(aSs) - 1
scalar p2 = rows(theta)
matrix theta_a = (p1>0) ? theta[1:p1]: {}
matrix A = mat_exp(theta_a, aSs)
matrix theta_b = (p2>p1) ? theta[p1+1:p2] : {}
matrix B = mat_exp(theta_b, bSs, 1)
</code>
</gretl-function>
<gretl-function name="PseudoHessAB" type="void" private="1">
<params count="4">
<param name="H" type="matrixref"/>
<param name="theta" type="matrixref"/>
<param name="dat" type="matrixref"/>
<param name="p1" type="int"/>
</params>
<code>matrix Sigma Ss A B
unscramble_dat(&dat, &Sigma, &Ss)
scalar p2 = rows(theta)
matrix aSs = Ss[,1:p1+1]
matrix bSs = Ss[,p1+2:p2+2]
scalar n = rows(Sigma)
scalar n2 = n*n
ABmat_exp(theta, aSs, bSs, &A, &B)
scalar ka = cols(aSs) - 1
scalar kb = cols(bSs) - 1
matrix S = zeros(2*n2, ka+kb)
if ka>0
S[1:n2,1:ka] = aSs[,1:ka]
endif
if kb>0
S[n2+1:,ka+1:ka+kb] = bSs[,1:kb]
endif
H = InfoMat(B, S, &A)
</code>
</gretl-function>
<gretl-function name="estAB" type="matrix" private="1">
<params count="9">
<param name="theta" type="matrixref"/>
<param name="Sigma" type="matrix"/>
<param name="aRd" type="matrix"/>
<param name="bRd" type="matrix"/>
<param name="vcv" type="matrixref" optional="true"/>
<param name="errcode" type="scalarref"/>
<param name="method" type="int"/>
<param name="verbose" type="int" default="1"/>
<param name="transfABkakb" type="matricesref" optional="true"/>
</params>
<code># new optional arg transfABkakb added by Sven for transfer/
# to avoid calc duplication
matrix aSs = imp2exp(aRd)
matrix bSs = imp2exp(bRd)
scalar npar = rows(theta)
scalar n = rows(Sigma)
scalar p1 = cols(aSs) - 1
scalar p2 = p1 + cols(bSs) - 1
matrix dat = vec(Sigma) ~ aSs ~ bSs
matrix g H
if verbose > 1
set max_verbose 1
else
set max_verbose 0
endif
scalar err = 0
if method == 0
catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, p1)")
err = $error
scoreAB(&g, &theta, &dat, p1)
elif method == 1
catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, p1)", "scoreAB(&g, &theta, &dat, p1)")
err = $error
elif method == 2
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, p1)")
err = $error
scoreAB(&g, &theta, &dat, p1)
elif method == 3
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, p1)", "scoreAB(&g, &theta, &dat, p1)")
err = $error
elif method == 4
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, p1)", "scoreAB(&g, &theta, &dat, p1)", "PseudoHessAB(&H, &theta, &dat, p1)" )
err = $error
endif
scalar crit = maxr(abs(g))
scalar warn = (crit > 1.0e-1)
if err || warn
matrix C = {}
errcode = err + 10*warn
outfile stderr --write
printf "err = %d; warn = %d; Gradient: %10.6f", err, warn, g
outfile --close
return C
endif
matrix A B
ABmat_exp(theta, aSs, bSs, &A, &B)
# use matlab-style "matrix division" for speed and accuracy
matrix C = A\B
if verbose > 1
printf "Estimated A matrix:\n%12.5f", A
printf "Estimated B matrix:\n%12.5f", B
printf "Estimated C matrix:\n%12.5f", C
printf "Check:\n%12.5f", cholesky(Sigma)
endif
# new transplanted block from SVAR_estimate():
n2 = n*n
ka = cols(aSs) - 1
kb = cols(bSs) - 1
matrix S = zeros(2*n2, ka+kb)
if ka > 0
S[1:n2,1:ka] = aSs[,1:ka]
endif
if kb > 0
S[n2+1:2*n2,ka+1:ka+kb] = bSs[,1:kb]
endif
if exists(vcv)
vcv = coeffVCV(S, &B, &A) # was coeffVCV(S, &mB, &mA)
endif
# transfer back the A, B, ka, kb results (ugly hack <Sven>)
if exists(transfABkakb)
transfABkakb = defarray(A, B, {ka}, {kb})
endif
return C
</code>
</gretl-function>
<gretl-function name="boot_printout" type="void" private="1">
<params count="5">
<param name="type" type="int"/>
<param name="n" type="int"/>
<param name="rep" type="int"/>
<param name="failed" type="int"/>
<param name="Spar_mat" type="matrixref"/>
</params>
<code>matrix bm = meanc(Spar_mat)
matrix bV = mcov(Spar_mat)
scalar nc = cols(Spar_mat)
scalar n2 = n*n
# force numerical zeros to real zeros
e = diag(bV) .< 1.0e-12
if maxc(e)
e = selifc(seq(1, nc), e')
bV[e,] = 0
bV[,e] = 0
endif
printf "Bootstrap results (%d replications, %d failed)\n", rep + failed, failed
if (type != 3) && ( cols(Spar_mat) == 2 * n2 )
# Long-run matrix at the end exists!
# And so this is a model with long-run restrictions (Bl-Quah-style or SVEC),
# or the user forced it via calc_lr.
matrix bK = mshape( bm[1:n2], n, n)
printStrMat(bK, bV[1:n2, 1:n2], "C")
matrix bL = mshape( bm[n2+1:], n, n)
printStrMat(bL, bV[1+n2: , 1+n2: ], "LongRun")
# endif
elif type != 3 # was: (type == 1) || (type == 2) || (type == 4)
# C model without printing the long-run matrix
# (SVEC / type 4 should not happen here, because it has long-run
# constraints by construction)
matrix bK = mshape(bm,n,n)
printStrMat(bK, bV, "C")
elif type == 3 # AB model
matrix bmA = mshape( bm[1:n2], n, n)
printStrMat(bmA, bV[1:n2, 1:n2], "A")
matrix bmB = mshape( bm[n2+1:], n, n)
printStrMat(bmB, bV[1+n2:,1+n2:], "B")
else
funcerr "shouldn't happen"
endif
</code>
</gretl-function>
<gretl-function name="max_eval" type="scalar" private="1">
<params count="1">
<param name="A" type="matrix"/>
</params>
<code>n = rows(A)
p = cols(A)/n
matrix compan = p==1 ? A : A | (I(n*(p-1)) ~ zeros(n*(p-1), n))
matrix lambda = eigengen(compan)
scalar maxmod = maxc(sumr(lambda.^2))
return maxmod
</code>
</gretl-function>
<gretl-function name="bias_correction" type="matrix" private="1">
<params count="7">
<param name="H" type="scalar"/>
<param name="Ahat" type="matrix"/>
<param name="mu" type="matrix" const="true"/>
<param name="E" type="matrix" const="true"/>
<param name="X" type="matrix" const="true"/>
<param name="Y0" type="matrix"/>
<param name="BC" type="matrixref"/>
</params>
<code>/* This function implements a bias correction for
the estimate of the VAR parameters as per Kilian, REStat (1998).
*/
n = rows(Ahat)
p = cols(Ahat)/n
k = cols(X)
cmu = cols(mu)
rY0 = rows(Y0)
#check for stationarity first
scalar maxmod = max_eval(Ahat)
if maxmod < 0.9999
matrix Ab = zeros(n, n*p)
loop i=1..H --quiet
matrix U = zeros(p,n) | resample(E)
if cmu > 0
U = mu + U
endif
matrix bY = varsimul(Ahat, U[rY0+1:,], Y0)
# printf "rows(bY) = %d, rows(X) = %d\n", rows(bY), rows(X)
matrix reg = X ~ mlag(bY, seq(1,p))
matrix Pi = mols(bY[p+1:,], reg[p+1:,])
matrix Ab += transp(Pi[k+1:k+n*p,])
endloop
Ab = Ab ./ H
matrix BC = Ahat - Ab
H = add_and_smash(&Ab, BC)
printf "H = %g\n", H
else
matrix Ab = Ahat
endif
return Ab
</code>
</gretl-function>
<gretl-function name="scoreC" type="void" private="1">
<params count="3">
<param name="ret" type="matrixref"/>
<param name="theta" type="matrixref"/>
<param name="dat" type="matrixref"/>
</params>
<code>matrix Sigma Ss
unscramble_dat(&dat, &Sigma, &Ss)
scalar n = rows(Sigma)
scalar npar = cols(Ss) - 1
matrix C = mat_exp(theta, Ss, 0)
matrix S = Ss[,1:npar]
matrix iC = inv(C)
matrix ret = iC' (qform(iC,Sigma) - I(n))
ret = vec(ret)'S
</code>
</gretl-function>
<gretl-function name="C1mat" type="matrix" private="1">
<params count="4">
<param name="A" type="matrix" const="true"/>
<param name="VECM" type="bool" default="0"/>
<param name="jalpha" type="matrix" optional="true" const="true"/>
<param name="jbeta" type="matrix" optional="true" const="true"/>
</params>
<code># (Jan 2018: used to be pointer args, coming from a time where it was necessary
# for optional args; now changed to non-pointers)
/*
computes C(1) out of the autoregressive matrices; jalpha and jbeta
are alpha and beta, which are used only if VECM is nonzero, that
is under cointegration
(Remark Sven: here jbeta was assumed to be (n,r)
without coeffs for restr. terms
Different from jbeta (n+1,r) in other contexts (?).
We now (Jan 2018) make sure only rows 1:n are used. */
scalar n = rows(A)
scalar p = cols(A) / n
matrix tmp = mshape(vec(A), n*n, p)
if VECM == 0
tmp = mshape(sumr(tmp), n, n)
matrix ret = inv(I(n) - tmp)
else # cointegrated
if !exists(jalpha) || !exists(jbeta)
funcerr "Need cointegration params for C1 in VECM!"
endif
matrix aperp = nullspace(jalpha')
matrix bperp = nullspace(jbeta[1:n, ]')
tmp = mshape(tmp*seq(1,p)', n, n)
matrix ret = bperp * inv(aperp'tmp*bperp) * aperp'
endif
return ret
</code>
</gretl-function>
<gretl-function name="init_C" type="matrix" private="1">
<params count="2">
<param name="Sigma" type="matrix"/>
<param name="Rd" type="matrix"/>
</params>
<code>matrix Ss = imp2exp(Rd)
scalar k = cols(Ss)
if k == 1
# nothing to estimate
ret = {}
else
scalar n = rows(Sigma)
matrix S = (k>1) ? Ss[,1:k-1] : {}
matrix s = Ss[,k]
matrix K = cholesky(Sigma)
matrix bigmat = (K ** I(n))
matrix ret = mols(vec(Sigma) - bigmat*s, bigmat*S)
endif
return ret
</code>
</gretl-function>
<gretl-function name="PseudoHessC" type="void" private="1">
<params count="3">
<param name="H" type="matrixref"/>
<param name="theta" type="matrixref"/>
<param name="dat" type="matrixref"/>
</params>
<code>matrix Sigma Ss
unscramble_dat(&dat, &Sigma, &Ss)
scalar n = rows(Sigma)
scalar npar = cols(Ss) - 1
matrix S = Ss[,1:npar]
# printf "PseudoHessC\n%10.4f\n%4.1f\n", theta, Ss
matrix C = mat_exp(theta, Ss, 0)
H = InfoMat(C, S) # was InfoMatC
</code>
</gretl-function>
<gretl-function name="estC" type="matrix" private="1">
<params count="7">
<param name="theta" type="matrixref"/>
<param name="Sigma" type="matrix"/>
<param name="Rd" type="matrix"/>
<param name="vcv" type="matrixref" optional="true"/>
<param name="errcode" type="scalarref"/>
<param name="method" type="int"/>
<param name="verbose" type="scalar" default="1"/>
</params>
<code>matrix Ss = imp2exp(Rd)
scalar n = rows(Sigma)
if cols(Ss) == 1
# nothing to estimate
matrix C = mshape(Ss, n, n)
if exists(vcv)
vcv = zeros(n*n,n*n)
endif
errcode = 0
return C
endif
scalar SCALE = 0 # EXPERIMENTAL
scalar npar = rows(theta)
if SCALE == 1
printf "Scale!\n"
matrix s = sqrt(diag(Sigma))
matrix sSig = Sigma ./ (s*s')
else
matrix sSig = Sigma
endif
# printf "estC\n%10.4f\n%4.1f\n", theta, Ss
matrix tmp = mat_exp(theta, Ss, 1)
if verbose > 2
/* obsolete ? */
printf "check within estC -- before estimation\n"
check_const(tmp, Rd)
endif
matrix dat = vec(sSig) ~ Ss
matrix g H
if verbose > 1
set max_verbose 1
else
set max_verbose 0
endif
err = 1
iters = 0
#set bfgs_toler 1.0e-03
matrix theta0 = theta
errcode = 0
loop while (err==1 && iters<100) --quiet
if method == 0
catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, -1)")
errcode = $error
scoreC(&g, &theta, &dat)
elif method == 1
catch scalar ll = BFGSmax(theta, "loglik(&theta, &dat, -1)", "scoreC(&g, &theta, &dat)")
errcode = $error
elif method == 2
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)")
errcode = $error
scoreC(&g, &theta, &dat)
elif method == 3
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)", "scoreC(&g, &theta, &dat)")
errcode = $error
elif method == 4
catch scalar ll = NRmax(theta, "loglik(&theta, &dat, -1)", "scoreC(&g, &theta, &dat)", "PseudoHessC(&H, &theta, &dat)")
errcode = $error
endif
if errcode>0
printf "errcode = %d\n", errcode
endif
scalar crit = maxr(abs(g))
err = (crit > 1.0e-4)
if err==1
iters++
theta = 0.1*mnormal(npar,1) + theta0
if verbose>1
printf "Iter %3d: Restarting... ll = %12.7f, crit = %16.10f\n", iters, ll, crit
printf "theta = %10.5f grad = %10.5f\n", theta', g
endif
endif
endloop
scalar crit = maxr(abs(g))
warn = (crit > 1.0e-1)
if !err
if (iters > 0) && (verbose > 1)
printf "Converged after %d restarts\n", iters
endif
# printf "estC(2)\n%10.4f\n%4.1f\n", theta, Ss
matrix C = mat_exp(theta, Ss, 1)
if SCALE == 1
C = C .* s'
endif
if verbose > 1
printf "Estimated C matrix:\n%12.5f", C
endif
if exists(vcv)
vcv = coeffVCV(Ss[,1:npar], &C)
endif
if verbose > 1
matrix Info = InfoMat(C, Ss[,1:npar]) # was InfoMatC
printf "estC : Info = \n%14.10f\n", Info
printf "estC : score = \n%14.10f\n", g
endif
else
if verbose > 1
printf "No convergence! :-( \n%12.6f\n", theta' | g
endif
matrix C = {}
errcode = 1
endif
return C
</code>
</gretl-function>
<gretl-function name="CheckNormalizeRd" type="scalar" private="1">
<params count="2">
<param name="R" type="matrixref"/>
<param name="d" type="matrixref"/>
</params>
<code>/*
Checks that
(1) the constraints are consistent
(2) the constraints are non-contradictory
if (1) fails, an error message is printed and R and d are replaced by
empty matrices; if (2) fails, redundant rows in R and d are dropped.
*/
p = rows(R)
r = rank(R)
ret = 0
if r < p
matrix Rd = R ~ d
if r < rank(Rd) #contradictory
R = {}
d = {}
ret = 1
else # redundant
matrix RR
matrix QQ = qrdecomp(Rd', &RR)
matrix e = abs(diag(RR)) .> $macheps
QQ = selifc(QQ, e')
matrix RR = selifr(selifc(RR, e'), e)
matrix Rd = QQ * RR
R = Rd[1:rows(Rd)-1,]'
d = Rd[rows(Rd),]'
ret = 2
endif
endif
return ret
</code>
</gretl-function>
<gretl-function name="add_constraint" type="scalar" private="1">
<params count="5">
<param name="Rd" type="matrixref"/>
<param name="n" type="int"/>
<param name="i" type="int"/>
<param name="j" type="int"/>
<param name="d" type="scalar"/>
</params>
<code>err = (i>n) || (j>n) || (i<1) || (j<1)
if !err
n2 = n*n
matrix tmp = zeros(1,n2 + 1)
k = n*(j-1) + i
tmp[1,k] = 1
tmp[1,n2+1] = d
# check for consistency/redundancy
if rows(Rd) > 0
matrix newR = Rd[,1:n2] | tmp[1:n2]
matrix newd = Rd[,n2+1] | d
else
matrix newR = tmp[1:n2]
matrix newd = d
endif
err2 = CheckNormalizeRd(&newR, &newd)
if err2==0
Rd = Rd | tmp
elif err2 == 1
printf "The restriction conflicts with previous ones and was ignored\n"
elif err2 == 2
printf "The restriction is redundant and was ignored\n"
endif
endif
return err ? -err : 10*err2 # -1 for bad input, 10 or 20 upstream error
</code>
</gretl-function>
<gretl-function name="cholRd" type="matrix" private="1">
<params count="1">
<param name="n" type="int"/>
</params>
<code>n2 = n*n
k = 1
matrix ret = {}
loop i=1..n -q
loop j=1..n -q
matrix tmp = zeros(1,n2+1)
if i > j
tmp[k] = 1
ret = ret | tmp
endif
k++
endloop
endloop
return ret
</code>
</gretl-function>
<gretl-function name="diag_Rd" type="matrix" private="1">
<params count="2">
<param name="n" type="int"/>
<param name="x" type="scalar"/>
</params>
<code>return selifr(I(n*n), vec(I(n))) ~ (x*ones(n,1))
</code>
</gretl-function>
<gretl-function name="free_diag_Rd" type="matrix" private="1">
<params count="1">
<param name="n" type="int"/>
</params>
<code>n2 = n*n
return selifr(I(n2)~zeros(n2,1), !vec(I(n)))
</code>
</gretl-function>
<gretl-function name="imp2exp" type="matrix" private="1">
<params count="1">
<param name="Rd" type="matrix"/>
</params>
<code>/*
Given the constraints in implicit form, returns the matrix [ S | s ]
of the constraints in explicit form
*/
p = cols(Rd)
matrix R = Rd[, 1:(p-1)]
matrix d = Rd[,p]
matrix s = R'invpd(R*R')*d
matrix S = nullspace(R | s')
return S ~ s
</code>
</gretl-function>
<gretl-function name="check_const" type="void" private="1">
<params count="2">
<param name="K" type="matrix"/>
<param name="fullRd" type="matrix"/>
</params>
<code>k = cols(fullRd) - 1
matrix R = fullRd[,1:k]
matrix d = fullRd[,k+1]
printf "Constraint matrix:\n"
print R
matrix tmp = R*vec(K)
print tmp
printf "Should be:\n"
print d
</code>
</gretl-function>
<gretl-function name="lrConstr" type="matrix" private="1">
<params count="2">
<param name="C1" type="matrix"/>
<param name="lrRd" type="matrix"/>
</params>
<code>if rows(lrRd) == 0
return {} # nothing to do
else
n = rows(C1)
k = cols(lrRd) - 1
matrix d = lrRd[,k+1]
return ( lrRd[,1:k] * (I(n) ** C1) ) ~ d
endif
return ret
</code>
</gretl-function>
<gretl-function name="get_full_Rd" type="matrix" private="1">
<params count="2">
<param name="obj" type="bundleref"/>
<param name="verbosity" type="int" default="0"/>
</params>
<code>scalar type = obj.type
scalar n = obj.n
matrix fullRd = obj.Rd1 # grab sr-restrictions first
# (for type 1 or type 2 with short-run only: nothing else to do)
lr_constr = rows(obj.Rd1l) ? 1 : 0 # further long-run constraints?
if type == 3
funcerr "This func not usable for AB models (yet)"
# because long-run restrictions aren't implemented for them
elif (type == 4) || lr_constr || obj.calc_lr
# generic C-model (includes SVEC) with some kind of long-run restr
# before ML estimation, we need to take into account the
# permanent/temporary shock classification in SVEC
if type == 4
matrix alpha = obj.jalpha # don't mix up with the scalar sig level obj.alpha!!
# trim beta from restr. exo terms if needed:
matrix beta = (obj.jcase % 2 == 0) ? obj.jbeta[1:n,] : obj.jbeta
r = cols(beta)
matrix C1 = C1mat(obj.VARpar, 1, alpha, beta)
matrix J = zeros(n-r, r) | I(r)
matrix ptRd = (J ** nullspace(alpha'))' ~ zeros(r*(n-r), 1)
fullRd = fullRd | ptRd # put the constraints together
else # can only be type 2 (or 1 if obj.calc_lr)
matrix C1 = C1mat(obj.VARpar, 0) # compute the lr matrix
endif
if lr_constr # avoid function call if unnecessary
fullRd = fullRd | lrConstr(C1, obj.Rd1l) # put the constraints together
endif
if verbosity > 1
matrix lrSigma = qform(C1, obj.Sigma) # compute the lr cov matrix
printf "Long-run matrix (C1): \n%8.3f\n", C1
printf "Long-run Sigma: \n%8.3f\n", lrSigma
endif
# store the C1 matrix for possible future use
matrix obj.C1 = C1
endif
return fullRd
</code>
</gretl-function>
<gretl-function name="SampleMatrices" type="matrices" private="1">
<params count="4">
<param name="Ra" type="matrix"/>
<param name="da" type="matrix"/>
<param name="Rb" type="matrix"/>
<param name="db" type="matrix"/>
</params>
<code>/*
Returns an array with A and B evaluated at a random point
in the parameter space (useful for checking the constraint
matrices).
Parameters: Ra, da, Rb, db = constraint matrices;
*/
n = sqrt(cols(Ra | Rb))
matrices ret
matrix tmp = imp2exp(Ra ~ da)
if cols(tmp) > 1
tmp *= muniform(cols(tmp)-1,1) | 1
endif
ret += mshape(tmp,n,n)
matrix tmp = imp2exp(Rb ~ db)
if cols(tmp) > 1
tmp *= muniform(cols(tmp)-1,1) | 1
endif
ret += mshape(tmp,n,n)
return ret
</code>
</gretl-function>
<gretl-function name="Dtn" type="matrix" private="1">
<params count="1">
<param name="n" type="int"/>
</params>
<code>/*
Creates the matrix \tilde{D}_n.
Output has n^2 rows and n*(n-1)/2 columns; any (n x n) skewsymmetric
matrix has a vectorised form which lies in the space spanned by the
columns of \tilde{D}_n
*/
p = round(n*(n-1)/2)
matrix A = zeros(1,n) | (lower(unvech(seq(1,p)')) ~ zeros(n-1,1))
matrix B = zeros(n^2,p)
matrix C
loop i=1..p --quiet
C = A.=i
B[,i] = vec(C .- C')
endloop
return B
</code>
</gretl-function>
<gretl-function name="Umat" type="matrix" private="1">
<params count="4">
<param name="n" type="int"/>
<param name="R" type="matrix"/>
<param name="d" type="matrix"/>
<param name="S" type="matrix"/>
</params>
<code># See Lucchetti(2006) -- left-hand side of matrix T;
# see eq. 26 on page 248
p = cols(S)
matrix ret = {}
loop i=1..p --quiet
ret |= R * (mshape(S[,i],n,n)' ** I(n))
endloop
return ret
</code>
</gretl-function>
<gretl-function name="Tmat" type="matrix" private="1">
<params count="4">
<param name="n" type="int"/>
<param name="R" type="matrix"/>
<param name="d" type="matrix"/>
<param name="S" type="matrix"/>
</params>
<code># See Lucchetti(2006) -- right-hand side of matrix T;
# see eq. 26 on page 248
p = cols(S)
matrix ret = {}
loop i=1..p --quiet
ret |= R * (I(n) ** mshape(S[,i],n,n))
endloop
return ret
</code>
</gretl-function>
<gretl-function name="strucond" type="scalar" private="1">
<params count="6">
<param name="n" type="int"/>
<param name="Ra" type="matrix"/>
<param name="da" type="matrix"/>
<param name="Rb" type="matrix"/>
<param name="db" type="matrix"/>
<param name="verbose" type="int" default="0"/>
</params>
<code>/*
Checks the structure condition and optionally
prints out some of the relevant matrices:
Parameters: Ra, da, Rb, db = constraint matrices;
iprint = output mode:
*/
matrix Sa = imp2exp(Ra~da)
matrix Sb = imp2exp(Rb~db)
if verbose > 1
print Ra da Rb db Sa Sb
endif
matrix Ua = Umat(n,Ra,da,Sa)
matrix Ub = Umat(n,Rb,db,Sb)
matrix Tb = Tmat(n,Rb,db,Sb)
if verbose > 1
print Ua Ub Tb
endif
matrix C = Ua ~ zeros(rows(Ua),(n*(n-1)/2))
C |= Ub ~ Tb * Dtn(n)
if verbose > 1
print C
endif
/* purge zero rows from C */
matrix e = maxr(abs(C)) .> 0
C = selifr(C,e)
if verbose > 1
printf "After filtering zero rows, C is %d x %d\n", rows(C), cols(C)
endif
matrix CC = C'C
if verbose > 1
print CC
endif
d = det(CC)
if d == 0
matrix nspace = nullspace(CC)
u_rank = cols(nspace)
if verbose
loop i=1..u_rank --quiet
printf "Q_%d = \n", i
printf "%6.1f", mshape(nspace[1:n*n,i],n,n)
printf "H_%d = \n", i
printf "%6.1f", mshape(Dtn(n) * nspace[n*n+1:,i],n,n)
endloop
endif
endif
return (d > 0)
</code>
</gretl-function>
<gretl-function name="ordercond" type="scalar" private="1">
<params count="4">
<param name="n" type="int"/>
<param name="Ra" type="matrix"/>
<param name="Rb" type="matrix"/>
<param name="verbose" type="int" default="0"/>
</params>
<code>/*
Checks the order condition and optionally
prints out some of the relevant matrices:
Parameters: Ra, Rb = constraint matrices;
iprint = output mode:
*/
p = 2*n*n - (n + 1)*n/2
if verbose > 0
printf "no. of constraints on A: %d\n", rows(Ra)
printf "no. of constraints on B: %d\n", rows(Rb)
printf "no. of total constraints: %d\n", rows(Ra) + rows(Rb)
printf "no. of necessary restrictions for the order condition = %d\n", p
endif
return rank(Ra)+rank(Rb) >= p
</code>
</gretl-function>
<gretl-function name="rankcond" type="scalar" private="1">
<params count="6">
<param name="n" type="int"/>
<param name="Ra" type="matrix"/>
<param name="da" type="matrix"/>
<param name="Rb" type="matrix"/>
<param name="db" type="matrix"/>
<param name="verbose" type="int" default="0"/>
</params>
<code>/*
Checks the rank condition as per Amisano-Giannini and optionally
prints out some of the relevant matrices:
Parameters: Ra, da, Rb db, = constraint matrices;
Note that in fact we check for the existence of solutions to
the system given in eq. (9), chapter 4. The condition discussed later
(matrix Q) is, sadly, wrong.
*/
matrices AB = SampleMatrices(Ra, da, Rb, db)
matrix A = AB[1]
matrix B = AB[2]
matrix BB = B*B'
matrix Q11 = Ra * (A' ** BB)
matrix Q21 = -Rb * (B' ** BB)
matrix Q22 = Q21 * Dtn(n)
matrix Q = (Q11 ~ zeros(rows(Ra), n*(n-1)/2)) | (Q21 ~ Q22)
scalar r = rank(Q)
if verbose > 1
loop foreach m Q11 Q21 Q22 Q --quiet
printf "\n$m:\n%7.2f", $m
endloop
endif
if verbose > 0
printf "rank condition: r = %d, cols(Q) = %d\n", r, cols(Q)
endif
return r == cols(Q)
</code>
</gretl-function>
<gretl-function name="ident" type="scalar" private="1">
<params count="5">
<param name="Ra" type="matrix"/>
<param name="da" type="matrix"/>
<param name="Rb" type="matrix"/>
<param name="db" type="matrix"/>
<param name="verbose" type="int" default="0"/>
</params>
<code>/*
Main function for checking identification.
The algorithm is described fully in Lucchetti (2006), "Identification Of Covariance Structures", Econometric
Theory, Cambridge University Press, vol. 22(02), p 235-257.
Parameters: Ra, da, Rb, db = constraint matrices.
*/
matrix locRa = Ra
matrix locda = da
matrix locRb = Rb
matrix locdb = db
# check on the constraints on A and B for inconsistencies
# and/or redundancies
err = CheckNormalizeRd(&locRa,&locda)
if err
printf "Contradictory constraints on A!\n"
return 0
endif
err = CheckNormalizeRd(&locRb,&locdb)
if err
printf "Contradictory constraints on B!\n"
return 0
endif
n = round(sqrt(cols(Ra | Rb)))
/* check for the order condition */
scalar id_o = ordercond(n, locRa, locRb, verbose)
string printout = id_o ? "OK" : "fails!"
printf "Order condition %s\n", printout
# /* check for the structure condition */
# scalar id_s = strucond(n, locRa, locda, locRb, locdb, verbose)
# string printout = id_s ? "OK" : "fails!"
# printf "Structure condition %s\n", printout
/* check for the rank condition */
scalar id_r = rankcond(n, locRa, locda, locRb, locdb, verbose)
string printout = id_r ? "OK" : "fails!"
printf "Rank condition %s\n", printout
return (id_o && id_r)
</code>
</gretl-function>
<gretl-function name="doIRF" type="void" private="1">
<params count="1">
<param name="SVARobj" type="bundleref"/>
</params>
<code>/*
constructs the structural VMA representation. Note
that the companion matrix is never used explicitely;
The output is not returned by the function, but rather
put into the bundle under the "IRFs" key.
*/
scalar type = SVARobj.type
matrix varA = SVARobj.VARpar
scalar H = SVARobj.horizon + 1
scalar n = SVARobj.n
if (type == 1) || (type == 2) || (type == 4)
matrix C = SVARobj.C # was: SVARobj.S1
elif type == 3
if inbundle(SVARobj, "C") # maybe not yet computed
matrix C = SVARobj.C
else
matrix C = SVARobj.S1 \ SVARobj.S2
endif
endif
matrix ret = zeros(H,n*n)
scalar np = SVARobj.p * n
matrix tmp = I(np)
matrix prd = zeros(np,np)
loop i=1..H --quiet
ret[i,] = vec(tmp[1:n,1:n] * C)'
if (np>n)
prd[n+1:np, ] = tmp[1:np-n, ]
endif
prd[1:n,] = varA * tmp
tmp = prd
endloop
if SVARobj.ncumul > 0
matrix to_cum = SVARobj.cumul
tmp = zeros(n,n)
tmp[to_cum,] = 1
sel = selifr(transp(seq(1,n*n)), vec(tmp))
ret[, sel] = cum(ret[, sel])
endif
SVARobj.IRFs = ret
</code>
</gretl-function>
<gretl-function name="check_bounds" type="scalar" private="1">
<params count="3">
<param name="s" type="int"/>
<param name="v" type="int"/>
<param name="n" type="int"/>
</params>
<code># negative s is allowed and means to flip the shock
if abs(s) > n || (s == 0)
return 1
elif (v > n) || (v < 1)
return 2
endif
return 0
</code>
</gretl-function>
<gretl-function name="IRFgrph" type="string" private="1">
<params count="10">
<param name="IRFmat" type="matrix"/>
<param name="snum" type="int"/>
<param name="vnum" type="int"/>
<param name="scale" type="scalar"/>
<param name="sname" type="string"/>
<param name="vname" type="string"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
<param name="cumul" type="matrixref" optional="true"/>
<param name="boot" type="bundleref" optional="true"/>
<param name="bc" type="int" min="0" max="2" default="0"/>
</params>
<code># bc: bias correction choice
scalar flip = (snum < 0)
scalar snum = abs(snum)
scalar n = round(sqrt(cols(IRFmat)))
if !isnull(boot)
scalar bootrep = boot.rep
scalar alpha = boot.alpha
else
scalar bootrep = 0
endif
# anything to cumulate?
scalar cumulate = !isnull(cumul) ? sumr(sumc(cumul .= vnum)) > 0 : 0
scalar tmpfname = int(1000*muniform(1,1))
string tmpfile = $windows ? sprintf("@dotdir\\irf%d.gp", tmpfname) : sprintf("@dotdir/irf%d.gp", tmpfname)
scalar k = (snum-1)*n + vnum
scalar h = rows(IRFmat)
matrix x = flip ? -IRFmat[,k]./scale : IRFmat[,k]./scale # matrix?
if bootrep > 0
matrix hicb = flip ? -boot.lo_cb : boot.hi_cb # matrix?
matrix locb = flip ? -boot.hi_cb : boot.lo_cb
matrix mdn = flip ? -boot.mdns : boot.mdns
matrix locb = locb[,k] ./ scale
matrix hicb = hicb[,k] ./ scale
matrix mdn = mdn[,k] ./ scale
scalar miny = minc((locb | x))
scalar maxy = maxc((hicb | x))
else
scalar miny = minc(x)
scalar maxy = maxc(x)
endif
scalar miny = (miny>0) ? 0 : miny
scalar maxy = (maxy<0) ? 0 : maxy
set force_decpoint on
outfile "@tmpfile" --write
printf "set yrange [%g:%g]\n", miny*1.05, maxy*1.05
printf "set xrange [%g:%g]\n", 0, (h-1)*1.05
printf "set xzeroaxis\n"
string l_sname = sname=="" ? sprintf("%d", snum) : sname
string l_vname = vname=="" ? sprintf("%d", vnum) : vname
printf "set title \"IRF: %s -> %s", l_sname, l_vname
if bc == 1
printf "; bias-correction = partial"
elif bc == 2
printf "; bias-correction = full"
endif
if cumulate > 0
printf " (cumulated)"
endif
printf "\"\n"
if bootrep > 0
if keypos == 0
printf "set key off\n"
elif keypos == 1
printf "set key outside\n"
elif keypos == 2
printf "set key below\n"
endif
printf "set style fill solid 0.125\n"
printf "plot '-' using 1:2:3 w filledcurve t 'Bstrap %d%% CI', \\\n", floor(100*alpha)
printf "'-' w l lw 1 lc 1 t 'Bstrap median', \\\n"
printf "'-' w l lw 2 lc -1 t 'IRF'\n"
loop i=1..h -q
printf "%d\t%g\t%g\n", i-1, locb[i,], hicb[i,]
endloop
printf "e\n"
loop i=1..h -q
printf "%d\t%g\n", i-1, mdn[i,]
endloop
printf "e\n"
else
printf "plot '-' w l lw 2\n"
endif
loop i=1..h -q
printf "%d\t%g\n", i-1, x[i]
endloop
printf "e\n"
outfile --close
set force_decpoint off
return tmpfile
</code>
</gretl-function>
<gretl-function name="FEVDgrph" type="string" private="1">
<params count="5">
<param name="FEVDmat" type="matrix"/>
<param name="v" type="scalar"/>
<param name="vname" type="string"/>
<param name="snames" type="strings"/>
<param name="keypos" type="int" min="0" max="2" default="1"/>
</params>
<code>n = round(sqrt(cols(FEVDmat)))
h = rows(FEVDmat) - 1
scalar tmpfname = int(10000*muniform(1,1))
if $windows
sprintf datfile "@dotdir\\fevd%d.txt", tmpfname
sprintf tmpfile "@dotdir\\fevd%d.gp", tmpfname
else
sprintf datfile "@dotdir/fevd%d.txt", tmpfname
sprintf tmpfile "@dotdir/fevd%d.gp", tmpfname
endif
matrix sel = (v-1)*n + seq(1,n)
set force_decpoint on
outfile "@tmpfile" --write
if keypos == 0
printf "set key off\n"
elif keypos == 1
printf "set key outside\n"
elif keypos == 2
printf "set key below\n"
endif
printf "set yrange [0:100]\n"
printf "set xrange [%g:%g]\n", 0, h
printf "set xzeroaxis\n"
printf "set title \"FEVD for %s\"\n", vname
printf "set style fill solid 0.25\n"
printf "set style histogram rowstacked\n"
printf "set style data histogram\n"
loop i=1..n --quiet
string sname = snames[i]
if i == 1
printf "plot '%s' using 2 t '%s', \\\n", datfile, sname
elif i == n
printf "\t'' using %d t '%s'\n", i+1, sname
else
printf "\t'' using %d t '%s', \\\n", i+1, sname
endif
endloop
outfile --close
outfile "@datfile" --write
printf "%12.4f\n", seq(0, h)' ~ 100*FEVDmat[,sel]
outfile --close
set force_decpoint off
return tmpfile
</code>
</gretl-function>
<gretl-function name="modelstring" type="string" private="1">
<params count="1">
<param name="type" type="int"/>
</params>
<code>string ret = ""
if type == 1
ret = "plain"
elif type == 2
ret = "C"
elif type == 3
ret = "AB"
elif type == 4
ret = "SVEC"
endif
return ret
</code>
</gretl-function>
<gretl-function name="modeltype" type="scalar" private="1">
<params count="1">
<param name="s" type="string"/>
</params>
<code>scalar ret = 0
string s2 = toupper(s) # case insensitive
if s2 == "PLAIN"
ret = 1
elif s2 == "C"
ret = 2
elif s2 == "AB"
ret = 3
elif s2 == "SVEC"
ret = 4
elif s2 == "KPSW"
print "Please use 'SVEC' for cointegrated SVARs,"
print "the code 'KPSW' is deprecated (but still works for now)."
ret = 4
endif
return ret
</code>
</gretl-function>
<gretl-function name="optimstring" type="string" private="1">
<params count="1">
<param name="type" type="int"/>
</params>
<code>string ret = ""
if type == 0
ret = "BFGS (numerical)"
elif type == 1
ret = "BFGS (analytical)"
elif type == 2
ret = "Newton-Raphson (numerical)"
elif type == 3
ret = "Newton-Raphson (analytical score)"
elif type == 4
ret = "Scoring algorithm"
endif
return ret
</code>
</gretl-function>
<gretl-function name="set_default_dimensions" type="matrix" private="1">
<params count="4">
<param name="mod" type="bundleref"/>
<param name="lY" type="list"/>
<param name="lX" type="list" optional="true"/>
<param name="varorder" type="int"/>
</params>
<code># trivial stuff: set up dimensions etc.
n = nelem(lY)
k = nelem(lX)
/* no of endogenous variables */
mod.n = n
/* no of exogenous variables */
mod.k = k
/* VAR order */
mod.p = varorder
/* horizon for IRFs etc */
mod.horizon = 10
if $pd == 4
/* quarterly: try 5 years */
mod.horizon = 20
elif $pd == 12
/* monthly: two years */
mod.horizon = 24
endif
/* sample size */
list everything = lY lX
smpl everything --no-missing
matrix mreg = { everything }
mod.T = rows(mreg)
mod.t1 = $t1
mod.t2 = $t2
return mreg
</code>
</gretl-function>
<sample-script filename="examples/simple_C.inp">
set verbose off
include SVAR.gfn
open sw_ch14.gdt
series infl = 400*ldiff(PUNEW)
rename LHUR unemp
list X = unemp infl
list Z = const
Mod = SVAR_setup("C", X, Z, 3)
Mod.horizon = 36
SVAR_restrict(&Mod, "C", 1, 2)
set stopwatch
SVAR_estimate(&Mod)
printf "Time (Cmodel) = %g\n", $stopwatch
fevdmat = FEVD(&Mod)
print fevdmat
#IRFplot(&Mod, 1, 1)
IRFsave("simple_C_11_noboot.pdf", &Mod, 1, 1)
set stopwatch
bfail = SVAR_boot(&Mod, 1024, 0.90)
printf "Failed = %d, Time (bootstrap) = %g\n", bfail, $stopwatch
IRFsave("simpleC.pdf", &Mod, 0, 0) # 0 means "all"
</sample-script>
</gretl-function-package>
</gretl-functions>
test.inp
Description: application/gretlscript
