proc (9)=bqvar(x,p,stps,dtrend,selector);
/*
Performs the Blanchard & Quah (1989) decompostion of a
unrestricted VAR and calculates the impulse-response
functions (IRFs), the variance decompostions (FEVDs) and
the historical decompositions (xr).

Syntax: {B,logl,IRFD,IRFL,FEVDD,FEVDL,e,u,xr} = BQVAR(x,p,stps,dtrend,selector);

Written by Martin W Johansson (date: Oct 31th, 2003)
                   National Institute of Economic Research
                   PO. BOX 3116
                   SE-103 62 Stockholm
                   martin_w.johansson AT konj DOT se

-------------------------------------------------------
The following input arguments are needed:

       x = A T*k matrix with the data, T=#obs, k=variables.
       p = Number of lags in unrestricted VAR.
    stps = Number of steps in IRFs and FEVDs.
  dtrend = 0,1,2 depending on how you want the determinstic
           components defined, see below...
       0 = No deterministic trends.
       1 = Intercept. (default)
       2 = Intercept and linear trend.
selector = A 1*k matrix where a zero (0) indicates a shock
that should be blanked out during HD and a one (1) 
indicates a variable that should be included. 
Ex. selector={1 0 1} would blank out the second shock. 
A final note: As the function assumes that the 
data is stationary you will usually input data in first
differnces and hence the HD will also be in first 
differences. To get the HD in levels you must accumulate.
-------------------------------------------------------
Output arguments: B,logl,IRFD,IRFL,FEVDD,FEVDL,e,u,xr

B     = matrix that transforms reduced form residuals to structural shocks.
logl  = log likelihood
IRFD  = IRFs in first diff.
IRFL  = IRFs in levels (i.e. accumulated).
FEVDD = FEVD in first diff.
FEVDL = FEVD in levels (i.e. accumulated).
e     = reduced form residuals. 
u     = structural shocks.
xr    = The historical decomposition. 


-------------------------------------------------------
NOTES: The function basically performs a lt Cholesky decompostion of the 
long-run VMA matrix. Since the factorisation is lt this means that the only
shock which may have long-run effects on the first variable is the first
shock while the first and second shock may have long-run effects on the second
variable. The first, second and third shock may have long-run effects on the 
third variable etc, etc...
The IRFs and FEVD are organised in 'variablewise'. Thus the first block of k 
columns holds the IRFs/FEVDs of the first variable to shock 1...k. The next block
of k columns holds the IRFs/FEVDs of the second variable to shock 1...k. 
Finally the orthogonalized (reduced form) residuals are stored in result.u (result.e)
-------------------------------------------------------
Reference: Blanchard & Quah (1989), 'The dynamic effects of
aggregate demand and supply disturbances', American Economic
Review, 79:4, 654-673.           
-------------------------------------------------------*/
local k, eyek, xorg, xlag, nobs, intercept, lintrend, a, at, ndf, const, trprmtr, e,
      sigma, logl, sumat, C1, invC1, G, B, u, J, AM, imp_resp, IRFD, IRFL, sum_fevd_D, 
      sum_fevd_L, FEVDD, FEVDL, er, xr;

k=cols(x);
eyek=eye(k);
xorg=x;

/*Estimating unrestrictred VAR*/
xlag={};
for i(1,p,1);
xlag=xlag~nclag(x,i,i);
endfor;
x=trimr(x,p,0);
xlag=trimr(xlag,p,0);
nobs=rows(xlag);
intercept=ones(nobs,1);
lintrend=seqa(1,1,nobs);

if dtrend==1;
xlag=intercept~xlag;
endif;
if dtrend==2;
xlag=intercept~lintrend~xlag;
endif;

a=x/xlag; 
at=a';
ndf=nobs-rows(at);

if dtrend==1;
at=a[2:rows(a),.]';
const=a[1,.]';
endif;
if dtrend==2;
at=a[3:rows(a),.]';
const=a[1,.]';
trprmtr=a[2,.]';
endif;

e=x-xlag*a;
sigma=(1/nobs)*e'*e;
logl=-0.5*nobs*(k*(1+ln(2*pi))+ln(det(sigma)));   /*Eviews 4.0 manual, p503*/
/*imposing long-run restrictions*/
if p==1;
sumat=at;
else;
sumat=reshape((sumc((reshape(vec(at),p,k*k)')')),k,k)';
endif;
C1=inv(eyek-sumat);		  
invC1=inv(C1);				     
G=chol(C1*sigma*C1')'; 						   
B=invC1*G;					     
u=(inv(B)*e')';

/*calculating the IRFs*/
J=eyek~zeros(k,k*p-k);
AM=(at|(eye((p-1)*k)~zeros((p-1)*k,k)));
imp_resp={};
for i(0,stps-1,1);
imp_resp=imp_resp|vec(((J*mexp(AM,i)*J')*B)')';
endfor;
IRFD=imp_resp;
IRFL=cumsumc(imp_resp);

/*calculating the FEVDs in differences*/
FEVDD=cumsumc(IRFD.^2);
sum_fevd_D=sumc(((ones(k,1).*.FEVDD).*(eyek.*.ones(stps,k)))');
FEVDD=100*FEVDD./((reshape(sum_fevd_D,k,stps)').*.ones(1,k));
/*calculating the FEVDs in Levels*/
FEVDL=cumsumc(IRFL.^2);
sum_fevd_L=sumc(((ones(k,1).*.FEVDL).*(eyek.*.ones(stps,k)))');
FEVDL=100*FEVDL./((reshape(sum_fevd_L,k,stps)').*.ones(1,k));

/*HD*/
er=(B*(selector.*u)')';
xr=xorg[1:p,.]|zeros(nobs,k); 
if dtrend==0;
for t(p+1,nobs+p,1);
xr[t,.]=(at*vec((rev((xr[t-p:t-1,.])))'))'+er[t-p,.];
endfor;
endif;

if dtrend==1;
for t(p+1,nobs+p,1);
xr[t,.]=const'+(at*vec((rev((xr[t-p:t-1,.])))'))'+er[t-p,.];
endfor;
endif;

if dtrend==2;
for t(p+1,nobs+p,1);
xr[t,.]=const'+trprmtr'+(at*vec((rev((xr[t-p:t-1,.])))'))'+er[t-p,.];
endfor;
endif;
retp(B,logl,IRFD,IRFL,FEVDD,FEVDL,e,u,xr);
endp;