Sven Schreiber 0.1 2011-06-30 Do permanent-transitory decomp in system, with uncertainty See Schreiber (2011, mimeo) endogenous lag order (levels) CI rank include constant? which perma-trans decomp? "Stock-Watson" "Gonzalo-Granger" bootstrap replications restricted exogenous unrestricted exogenous /* this function determines the confidence bands for all variables endo: list of endogenous variables in Vecm rexo: list of restricted exog. vars (should be demeaned...) uexo: list of unrestricted exog. vars (no demeaning needed?) p: lag order r: cointegration rank withconst: include unrestricted constant in Vecm? coverage: 1-alpha error of the confidence band numrep: number of bootstrap replications PTtype: "SW" or "GG" */ string PTtype = (PTtypecode=0) ? "SW" : "GG" error = errorcheck(endo,r,PTtype,0) n = nelem(endo) numrexo = nelem(rexo) numuexo = nelem(uexo) matrix endomat = {endo} matrix rexomat = {rexo} matrix uexomat = {uexo} # save initial values matrix initvals = endomat[1:p,] # pxn smpl +p +0 # to make sure we never reach back to pre-sample starting values # estimate and save stuff from the original model matrix betaext = {} matrix shortruncoeffs = {} matrix Omega = {} matrix vecmresmat = {} matrix companmat = {} Teff = 0 estimateVecm(endo,rexo,uexo,p,r,withconst, &betaext,&shortruncoeffs,&Omega, &vecmresmat,&companmat,&Teff) matrix alphaT = shortruncoeffs[n*(p-1)+withconst+numuexo+1:,] #rxn if withconst matrix muT = shortruncoeffs[1,] # const term 1xn endif if numuexo matrix muuexoT = shortruncoeffs[n*(p-1)+1+withconst: n*(p-1)+withconst+numuexo,] endif if numrexo matrix betaextrexo = betaext[n+1:,] endif # construct the trans component (point estimate) matrix ycondmat = getycondmat(endo,rexo,uexo,p,withconst,PTtype) matrix GorS = getGorS(shortruncoeffs,betaext,p,withconst,PTtype) matrix gapmat = GorS*ycondmat # nxT # prepare the result matrices loop i=1..n --quiet matrix gaptilde$i = {} # will be TxW endloop loop numrep # bootstrap replications # 1) resample the innovations matrix newresmat = resample(vecmresmat) # Txn? # 2) construct artificial data # 2a combine new shocks, constant, and unrestr. exog. matrix tempmat = newresmat if withconst matrix tempmat += ones(Teff,1)*muT endif if numuexo matrix tempmat += uexomat[p+1:,]*muuexoT endif if numrexo matrix tempmat += rexomat[p+1:,]*betaextrexo*alphaT endif # 2b simulate matrix newdatamat = varsimul(companmat[1:n,],tempmat,initvals) # 3) estimate on the new data list newdata = null smpl -p +0 # restore the original sample to create series including init vals loop i=1..n --quiet series temp$i = newdatamat[,i] # ($i because names in list must be different) list newdata += temp$i endloop smpl +p +0 # and again avoid reach-back to pre-sample values if withconst vecm p r newdata; uexo; rexo --silent else vecm p r newdata; uexo; rexo --nc --silent endif matrix newbeta = $jbeta matrix newshruncoeffs = $coeff # 4) construct the trans component matrix newGorS = getGorS(newshruncoeffs,newbeta,p,withconst,PTtype) matrix newgapmat = newGorS*ycondmat # nxT # 5) store the results (of all variables) loop i=1..n --quiet matrix gaptilde$i ~= transp(newgapmat[i,]) endloop endloop # end bootstrap loop # evaluate the results list out = null loop i=1..n --quiet # A) naive bootstrap: series @PTtype_naivelow$i = (quantile((gaptilde$i)',(1-coverage)/2))' series @PTtype_naiveup$i = (quantile((gaptilde$i)',1-(1-coverage)/2))' # B) and Hall's band: # subtract the point estimate: matrix tilHall$i = gaptilde$i- (gapmat[i,])' **ones(1,numrep)) # quantile of this "centered" dist: matrix Htemplow$i = (quantile((tilHall$i)',(1-coverage)/2))' matrix Htempup$i = (quantile((tilHall$i)',1-(1-coverage)/2))' # and "switch" it around: series @PTtype_Halllow$i = (gapmat[i,])' - Htempup$i series @PTtype_Hallup$i = (gapmat[i,])' - Htemplow$i # construct the output series @PTtype_gap$i = gapmat[i,]' list out += @PTtype_gap$i @PTtype_naivelow$i @PTtype_naiveup$i @PTtype_Halllow$i @PTtype_Hallup$i endloop return out endogenous lag order (levels) CI rank include constant? which perma-trans decomp? "Stock-Watson" "Gonzalo-Granger" apply to which variable (number in list) restricted exogenous unrestricted exogenous /* as indicated by 'whichvar', this function only determines the confidence bands for a single variable endo: list of endogenous variables in Vecm rexo: list of restricted exog. vars (should be demeaned...) uexo: list of unrestricted exog. vars (no demeaning needed?) p: lag order r: cointegration rank withconst: include unrestricted constant in Vecm? whichvar: choose variable number for which to calculate stuff coverage: 1-alpha error of the confidence band SWP: flag whether SWP decomp is used (else GG) */ string PTtype = (PTtypecode=0) ? "SW" : "GG" error = errorcheck(endo,r,PTtype,whichvar) numrexo = nelem(rexo) matrix betaext = {} matrix shortruncoeffs = {} matrix Omega = {} estimateVecm(endo,rexo,uexo,p,r,withconst, &betaext,&shortruncoeffs,&Omega) matrix vecshrun = vec(shortruncoeffs) # point estimate of gap matrix irow = get_rowi(vecshrun,betaext,p,numrexo,withconst,whichvar,PTtype) # the Jacobian for the i-th row: matrix Ji = fdjac(vecshrun, get_rowi(vecshrun,betaext,p,numrexo,withconst,whichvar,PTtype)) matrix ycondmat = getycondmat(endo,rexo,uexo,p,withconst,PTtype) matrix gapimat = irow'ycondmat # standard error of the gap matrix gapisemat = sqrt(diag(qform(ycondmat'Ji,Omega))) # transfer to a series output sprintf i "%d", whichvar # to get the string rep of whichvar se_multiplier = invcdf(N,(1-coverage)/2) # negative series @PTtypegap@i = gapimat series @PTtype@i_Delt_se = gapisemat series @PTtype@i_Delt_low = @PTtypegap@i+se_multiplier*@PTtype@i_Delt_se series @PTtype@i_Delt_hi = @PTtypegap@i-se_multiplier*@PTtype@i_Delt_se list out = @PTtypegap@i @PTtype@i_Delt_se @PTtype@i_Delt_low @PTtype@i_Delt_hi return out if withconst vecm p r endo; uexo; rexo --silent else vecm p r endo; uexo; rexo --nc --silent endif matrix betaext = $jbeta matrix shortruncoeffs = $coeff # full cov. of vec(shortruncoeffs): matrix Omega = $sigma**$xtxinv # extended stuff (only needed for bootstrap version): if !isnull(Teff) Teff = $T endif if !isnull(vecmresmat) matrix vecmresmat = $uhat # Txn endif if !isnull(companmat) matrix companmat = $compan endif /* Build the matrix that gives the SWP trans component from gretl's VECM output. This must match with getycondmat(). The dimensions of the return matrix should depend on the lag order and whether an unrestricted constant is included or not. shortruncoeffs: $coeff from gretl's Vecm estimate betaext: $jbeta, including coeffs. from restr. exo vars p: lag order (levels) withconst: flag if (unrestricted) constant is included The gretl ordering of the coeffs is: 1st row: constant 2nd row: first var, first diff-lag 3rd row: first var, second diff-lag ... p th row: first var, (p-1)th diff-lag ... (i-1)*(p-1)+1+j th row: i-th var, j-th diff lag ... n*(p-1)+1+withconstant th row: first unrestr. exo ... n*(p-1)+withconstant+numuexo th row: last unrestr. exo n*(p-1)+1+withconstant+numuexo th row: first ECT loadings ... n*(p-1)+withconstant+numuexo+r th row: last ECT loadings */ n = cols(shortruncoeffs) # how many endogenous variables r = cols(betaext) # cointegration rank numuexo = rows(shortruncoeffs)-n*(p-1)-withconst-r matrix beta = betaext[1:n,] # only coeffs of endo vars # alphas (could also be retrieved by $jalpha...) matrix alpha = transp(shortruncoeffs[n*(p-1)+1+withconst+numuexo:,]) /* here we don't use getB1() because we also need the side product of the B$j matrices here */ matrix Bsum = zeros(n,n) # all B-matrices added loop j=1..p-1 --quiet # which lag # very clever indexing: matrix idx = seq(withconst+j,(n-1)*(p-1)+withconst+j, p-1) matrix B$j = transp(shortruncoeffs[idx,]) matrix Bsum += B$j # useful for B(1) endloop matrix B1 = I(n)-Bsum matrix B1abI = inv(B1-alpha*beta') # (\Gamma(1)-\alpha\beta')^{-1} matrix b_aI = inv(beta'B1abI*alpha) # just an abbrev. matrix P = B1abI*alpha*b_aI*beta' # from Hecq et al. matrix Pext = B1abI*alpha*b_aI*betaext' # with restr. exo matrix GorS = Pext # starting point if PTtype="SW" # the B* matrices (and also their trafo) matrix Bstarsum = zeros(n,n) loop j=0..p-2 --quiet matrix Bstar$j = zeros(n,n) loop i=j+1..p-1 --quiet matrix Bstar$j += B$i endloop matrix Bstarsum += Bstar$j endloop loop j=0..p-2 --quiet matrix GorS ~= -(I(n)-P)*B1abI*Bstar$j endloop endif # now put together the other parts ## smu / constant term / unrestricted exo vars if withconst || numuexo matrix gorsmu = B1abI*alpha*b_aI*b_aI*beta' if PTtype="SW" matrix gorsmu += (I(n)-P)*B1abI*Bstarsum*(I(n)-P) endif # new with out-factored muext: matrix muext = get_muext(shortruncoeffs,p,withconst,numuexo) matrix gorsmuext = gorsmu*B1abI*muext matrix GorS ~= gorsmuext # else S stays the same, not extended endif return GorS /* construct the y_cond vector (rather the corresponding data matrix) for the GG/SWP trans comp; must match getG()/getS() */ numuexo = nelem(uexo) matrix endomat = {endo} matrix rexomat = {rexo} matrix uexomat = {uexo} matrix ycondmat = endomat~rexomat if PTtype="SW" # add lagged differences for SW matrix diffp2lag = mlag(diff(endomat),seq(0,p-2)) matrix ycondmat = ycondmat~diffp2lag endif # constant term and unrestricted exo vars matrix ycondmat = withconst ? ycondmat~ones(rows(endomat),1) : ycondmat matrix ycondmat = numuexo ? ycondmat~uexomat : ycondmat matrix ycondmat = transp(ycondmat) return ycondmat error = 0 # check valid PTtype input if PTtype="SW" || PTtype="GG" # ok, pass else error = 1 funcerr "PTtype must be SW or GG" endif # check validity of CI rank n = nelem(endo) if r>n-1 error = 1 funcerr "cointegration rank exceeds n-1" endif # check bounds of whichvar if whichvar>n error = 1 funcerr "whichvar out of bounds" endif # rule out missing data series miss = missing(endo) if sum(miss)>0 error = 1 funcerr "found missing data, cannot handle that" endif return error /* construct the constant coeff and unrestricted exo coeffs -- for the structure of shortruncoeffs, see getS() */ wc = withconst # just an abbrev n = cols(shortruncoeffs) # how many endogenous variables # constant term (still comes first) matrix mu = wc ? transp(shortruncoeffs[1,]) : {} # other unrestricted terms matrix muuexo = numuexo ? transp(shortruncoeffs[n*(p-1)+1+wc:n*(p-1)+wc+numuexo,]) : {} # put it together matrix muext = mu~muuexo return muext /* transform the VECM coefficients that are spit out by gretl to the (transpose of the) row of G or S (g'_i or s'_i) that will give the trans SWP/GG component. This function is the one for which we need the Jacobian. vecshrun: vec($coeff) from gretl's Vecm estimate beta: $jbeta, cointegration coefficients p: lag order (levels) numrexo: number of restricted exo withconst: whether an unrestricted constant is included whichvar: indicator to pick which row of S PTtype: "SW" or "GG" */ n = rows(betaext)-numrexo numrows = rows(vecshrun)/n # rows of coeffmat matrix shortruncoeffs = mshape(vecshrun,numrows,n) # unvec matrix GorS = getGorS(shortruncoeffs,betaext,p,withconst,PTtype) return transp(GorS[whichvar,]) testdelta = 1 # apply the delta method testboot = 1 # and/or apply the bootstrap ### test code, function call at end of file: open AWM.gdt --preserve # YER # real euro GDP # ITR # real euro gross investment # PCR # real private consumption string type = "SW" # Stock-Watson "SW" or Gonzalo-Granger "GG" list endo = YER ITR PCR whichvar = 1 # only relevant for delta method p = 4 r = 2 withconst = 1 numrep = 500 coverage = 0.9 if !islist(rexo) list rexo = null endif if !islist(uexo) list uexo = null endif if testdelta list deltagapstuff = getDeltaout(endo,p,r,withconst,type,\ coverage,whichvar,rexo,uexo) endif if testboot list bootgapstuff = getBootout(endo,p,r,withconst, \ type,coverage,numrep,rexo,uexo) endif