laubach and williams code

The error "matrix not positive definite" happens when you are trying to invert a matrix using the function, invpd. The invpd function requires that the input matrix is positive definite. Most of the time it is used to invert a moment matrix (X'X) which should be positive definite unless there is some linear dependence in the variables.

When you say your matrix is full of zeros, do you mean that it is entirely zeros? Or that it just has quite a few of them? You might be having a problem loading your data correctly. Can you post a link to the code so we can offer more help?

Thank you for answering me.

I mean when i use the function "print" to see my matrix. I see that it's entirely zeros

If your matrix is entirely zeros, then it is certainly singular. Does the program simulate data, or is this data that you loaded from a file?

thank you again

well,

here is my code

/* ----------------------------------------------------------------------
rstar_s1.g: First-stage estimation of the model without hours in
Laubach-Williams, FEDS 2001-56. Last modified 09/03/02.
---------------------------------------------------------------------- */

library optmum;
optset;
output file = /mq/home/m1hxl00/work/nri/nrv2/n2v1.out reset;
format /rd 8,4;

load dat[218,7] = /mq/home/m1hxl00/work/nri/nrv2/nr2_0902.dat;
dat[.,5] = dat[.,5]-dat[.,3]; /* Relative oil price inflation */
dat[.,6] = dat[.,6]-dat[.,3]; /* Relative import price inflation */
dat[.,7] = dat[.,7]-dat[.,4]; /* Ex ante real rate */
mv = miss(0,0);

nt = 166; /* Sample 61:1 - 02:2 */
nar = 7; /* Number of regressors */
nsv = 3; /* Number of state variables */
nst = nt; /* Startup periods for filter */
niter = 1000; /* Number of draws for generating error bands */
seed1 = 40;
indx = 0; /* Index for initialisation of P(0|0) during MC simulation */

/* Setting up system matrices:
a1 filtered state a(t|t)
a2 smoothed state a(t|T)
p1s covariance matrix of filtered state
p2s covariance matrix of smoothed state
b coefficients on explanatory variables in measurement equation
xr regressor matrix in measurement equation
z coefficients on state vector in measurement equation
h covariance matrix of measurement equation error
t transition matrix in transition equation
q covariance matrix of transition equation error */

a1 = zeros(nt+1,nsv);
a2 = zeros(nt+1,nsv);
p1s = zeros(nst*nsv,nsv);
p2s = zeros(nst*nsv,nsv);
b = zeros(2,nar);
z = zeros(2,nsv);
h = zeros(2,2);
t = zeros(nsv,nsv);
t[1:2,1] = ones(2,1);
t[3,2] = 1;
q = zeros(nsv,nsv);

/* Computing starting values for parameters */

/* Potential XBO proxied by segmented trend 60:1 to 00:4 */
xr = ones(nt+4,1)~seqa(1,1,nt+4);
xr = xr~(zeros(56,1)|seqa(1,1,nt-52))~(zeros(142,1)|seqa(1,1,nt-138));
xi = invpd(xr'xr);
bpot = xi*xr'dat[49:nt+52,1];
gpot = xr*bpot;
gap = 100*(dat[49:nt+52,1]-xr*bpot);

/* IS equation */
yis = gap[5:nt+4];
xis = gap[4:nt+3]~gap[3:nt+2];
//******************************************
//Line below fails
bis = invpd(xis'xis)*xis'yis;
ris = yis-xis*bis;
sis = sqrt(ris'ris/(nt-cols(xis)));

/* Phillips equation */
yph = dat[53:nt+52,3];
xph = dat[52:nt+51,3]~(dat[51:nt+50,3]+dat[50:nt+49,3]+dat[49:nt+48,3])/3;
xph = xph~(dat[48:nt+47,3]+dat[47:nt+46,3]+dat[46:nt+45,3]+dat[45:nt+44,3])/4;
xph = xph~gap[4:nt+3]~dat[52:nt+51,5]~dat[53:nt+52,6];
//******************************************
//Line below fails
bph = invpd(xph'xph)*xph'yph;
rph = yph-xph*bph;
sph = sqrt(rph'rph/(nt-cols(xph)));

/* Transition equation std devs */
s3 = 0.5;
d = 0.85;

/* Constructing data matrices and vector of initial values */
y = 100*dat[53:nt+52,1]~yph;
xr = 100*(dat[52:nt+51,1]~dat[51:nt+50,1])~xph[.,1:3]~xph[.,5:6];

startv = bis|bph[1:2]|bph[4:6]|d|sis|sph|s3;

output on;
"Starting values";
startv';
"";
output off;

err = reshape(mv,nt,2);
lf = reshape(mv,nt,1);
{val,f,g,retcode} = optmum(&li,startv);

/* res stores smoothed state, ref stores filtered state */

res = a2[2:nt+1,1];
ref = a1[2:nt+1,1];

/* Initial smoothed state and its covariance matrix are stored for
generating random draws during simulations */
stin = a2[2,.];
pp1 = p2s[1:nsv,.];
{tx,ty} = eigrs2(pp1);
dg = diagrv(eye(rows(tx)),tx);
hh2 = ty*sqrt(dg);

/* Computing gradients */
llf = lf;
dval = maxc(10^(-6)*abs(val')|ones(1,rows(val))*10^(-6));
bval = reshape(val,rows(val),rows(val))';
bval = diagrv(bval,val+dval);
i = 1;
do while i <= rows(val);
    call li(bval[.,i]);
    llf = llf~lf;
    i = i+1;
endo;
llf = (llf[.,2:rows(val)+1]-llf[.,1])./dval';
info = invpd(llf'llf);
bse = sqrt(diag(info));

output on;
"Likelihood value";
-f;
"";
"MLE Parameters of Model, SE, t Statistics";
val~bse~(abs(val)./bse);
"";
"Sum of coeffs on gap in IS equation";
sumc(val[1:2]);
"";

"Box-Pierce Q Statistics, L = 1,...,16, and 5% critical values";
bpq = zeros(16,2);
i = 1;
do while i <= 16;
    bpq[i,1] = ((err[1+i:nt,1]'err[1:nt-i,1])/(err[.,1]'err[.,1]))^2;
    bpq[i,2] = ((err[1+i:nt,2]'err[1:nt-i,2])/(err[.,2]'err[.,2]))^2;
    i = i+1;
endo;

i = 2;
do while i <= 16;
    bpq[i,.] = bpq[i,.]+bpq[i-1,.];
    i = i+1;
endo;
cv = {3.84,5.99,7.82,9.49,11.07,12.59,14.07,15.51,16.92,18.31,
    19.68,21.03,22.36,23.69,25,26.3};
nt*bpq~cv;
"";

/* Computing error bands for state vector using Monte Carlo */

indx = 1; /* Switches computation of initial state from GLS to draw */
cum1 = zeros(nt,1);
cum2 = zeros(nt,1);
{tx,ty} = eigrs2(info);
dg = diagrv(eye(rows(tx)),tx);
hh = ty*sqrt(dg);
i = 1;
do while i <= niter;
    theta = hh*rndns(rows(val),1,seed1)+val;
    call li(theta);
    cum1 = cum1+(a2[2:nt+1,1]-res[.,1]).*(a2[2:nt+1,1]-res[.,1]);
    j = 1;
    do while j <= nst-1;
        cum2[j] = cum2[j]+p2s[j*nsv+1,1];
        j = j+1;
    endo;
    cum2[nst:nt] = cum2[nst:nt]+ones(nt-nst+1,1)*p1s[(nst-1)*nsv+1,1];
    i = i+1;
endo;
cum1 = cum1/niter;
cum2 = cum2/niter;
bnd2 = sqrt(cum1+cum2);

"Mean SE";
meanc(bnd2);
(100*dat[53:nt+52,1]-res[.,1])~bnd2;
res/100;

end;
/* ----- Define the log likelihood function for each observation ------*/
proc li(inp);
    local p1,p2,i,f,fi,k,ldf,vv,j,flag1,c,x,om,yy,tmp;
    
    /* Insert parameters into system matrices */
    
    a1 = zeros(nt+1,nsv);
    a2 = zeros(nt+1,nsv);
    p1s = zeros(nst*nsv,nsv);
    p2s = zeros(nst*nsv,nsv);
    b = zeros(2,cols(xr));
    b[1,1:2] = inp[1:2]';
    b[2,1] = inp[5];
    b[2,3:4] = inp[3:4]';
    b[2,5] = 1-sumc(inp[3:4]);
    b[2,6:7] = inp[6:7]';
    z[1,1] = 1;
    z[1,2:3] = -inp[1:2]';
    z[2,2] = -inp[5];
    c = inp[8]|zeros(nsv-1,1);
    h[1,1] = inp[9]^2;
    h[2,2] = inp[10]^2;
    q[1,1] = inp[11]^2;
    
    /* Starting values for a(0|0), P(0|0) */
    
    if indx eq 1;
        p1 = pp1;
        a1[1,.] = (hh2*rndns(nsv,1,seed1)+stin')';
        goto flag1;
    endif;
    
    x = z|z*t|z*t*t|z*t*t*t;
    om = zeros(8,8);
    
    /* First initialize lower triangular off-diagonal elements */
    om[5:6,3:4] = z*t*q*z';
    om[7:8,3:4] = z*t*t*q*z';
    om[7:8,5:6] = z*(t*t*q*t'+t*q)*z';
    /* Initialize upper triangular off-diagonal elements by adding transpose */
    om = om+om';
    /* Initialize diagonal elements */
    om[1:2,1:2] = h;
    om[3:4,3:4] = z*q*z'+h;
    om[5:6,5:6] = z*(t*q*t'+q)*z'+h;
    om[7:8,7:8] = z*(t*t*q*t'*t'+t*q*t'+q)*z'+h;
    
    if det(om) le 10^(-6);
        "om" det(om);
        a1[1,.] = 100*gpot[5:3]';
        p1 = 2*eye(nsv);
        goto flag1;
    endif;
    om = invpd(om);
    if det(x'om*x) le 10^(-6);
        "x'om*x" det(x'om*x);
        a1[1,.] = 100*gpot[5:3]';
        p1 = 2*eye(nsv);
        goto flag1;
    endif;
    p1 = invpd(x'om*x);
    yy = y[1,.]'|y[2,.]'|y[3,.]'|y[4,.]';
    tmp = b*xr[1,.]'|(b*xr[2,.]'+z*c)|(b*xr[3,.]'+z*c+z*t*c);
    tmp = tmp|(b*xr[4,.]'+z*(eye(nsv)+t+t*t)*c);
    a1[1,.] = (p1*x'om*(yy-tmp))';
    tmp = yy-tmp-x*a1[1,.]';
    p1 = p1*(tmp'tmp)/3;
    
    flag1:
        
        lf = zeros(nt,1);
    
    /* Starting the Kalman filter:
    a1 stores a(t|t), t = 0,...,nt, a2 stores a(t|t-1), t = 1,...,nt
    p1s stores P(t|t), p2s stores P(t|t-1), t = 1,...,nst */
    
    i = 1;
    do while i <= nst;
        a2[i,.] = a1[i,.]*t'+c';
        err[i,.] = y[i,.]-a2[i,.]*z'-xr[i,.]*b';
        p2 = t*p1*t'+q;
        p2s[(i-1)*nsv+1:i*nsv,.] = p2;
        f = z*p2*z'+h;
        fi = inv(f);
        k = p2*z'fi;
        a1[i+1,.] = a2[i,.]+err[i,.]*k';
        p1 = (eye(nsv)-k*z)*p2;
        p1s[(i-1)*nsv+1:i*nsv,.] = p1;
        lf[i,1] = -0.5*(cols(y)*ln(2*pi)+ln(det(f))+err[i,.]*fi*err[i,.]');
        err[i,.] = err[i,.]./sqrt(diag(f))'; /* end filter */
        i = i+1;
    endo;
    
    ldf = ln(det(f));
    
    /* Begin steady state kalman filter */
    i = nst+1;
    do while i <= nt;
        a2[i,.] = a1[i,.]*t'+c';
        err[i,.] = y[i,.]-a2[i,.]*z'-xr[i,.]*b';
        a1[i+1,.] = a2[i,.]+err[i,.]*k';
        lf[i,1] = -0.5*(cols(y)*ln(2*pi)+ln(det(f))+err[i,.]*fi*err[i,.]');
        err[i,.] = err[i,.]./sqrt(diag(f))'; /* end filter */
        i = i+1;
    endo;
    
    /* Smoothing
    In a2, a(t|t-1) is overwritten by a(t-1|T)
    In p2s, P(t|t-1) is overwritten by P(t-1|T) */
    
    /* Recursively with values from steady state filter for t = nt,...,nst */
    
    j = p1s[(nst-2)*nsv+1:(nst-1)*nsv,.]*t'inv(p2s[(nst-1)*nsv+1:nst*nsv,.]);
    a2[nt+1,.] = a1[nt+1,.];
    
    i = 1;
    do while i <= nt+1-nst;
        a2[nt+1-i,.] = (a1[nt+1-i,.]'+j*(a2[nt+2-i,.]'-a2[nt+1-i,.]'))';
        i = i+1;
    endo;
    
    /* Recursively with values from startup period for t = nst-1,...,1 */
    
    p2s[(nst-1)*nsv+1:nst*nsv,.] = p1s[(nst-2)*nsv+1:(nst-1)*nsv,.]
        +j*(p1s[(nst-1)*nsv+1:nst*nsv,.]-p2s[(nst-1)*nsv+1:nst*nsv,.])*j';
    i = 1;
    do while i <= nst-2;
        j = p1s[(nst-2-i)*nsv+1:(nst-1-i)*nsv,.]
            *t'inv(p2s[(nst-1-i)*nsv+1:(nst-i)*nsv,.]);
        a2[nst-i,.] = (a1[nst-i,.]'+j*(a2[nst+1-i,.]'-a2[nst-i,.]'))';
        p2s[(nst-1-i)*nsv+1:(nst-i)*nsv,.] = p1s[(nst-2-i)*nsv+1:(nst-1-i)*nsv,.]
            +j*(p2s[(nst-i)*nsv+1:(nst+1-i)*nsv,.]
            -p2s[(nst-1-i)*nsv+1:(nst-i)*nsv,.])*j';
        i = i+1;
    endo;
    
    vv= sumc(lf);
    retp(-vv);
    
endp;

it's not working in the lines that are marked as failing. Basically i don't understand why this matrix is entirely zero

thank again

What matrix is entirely zero? Do you mean xph and xis, or the matrices made from xph'xph, xis'xis, or another one? Is the data set you are using an example dataset that came with the code, or is it your data?

I mean in the IS equation, xis and yis, when i print i just see 2 columns of zero

Normally it's not suppose to do

i use data coming from a excel sheet

thanks so much for taking time to answer me

I know it is not the problem, but you don't want to invert a matrix because it's numerically not a good idea. I recommend the following. Instead of inv(X) or invpd(X), you should write eye(n)/X. This avoids the numerical problems of inverting a matrix.