trouble with code2

This code below was posted by another user. But I am using the same code generally speaking.

At the bottom of the code where the long results are computed appear to be exactle the same as the short run.

I have tried to figure out why this happens but can't find the source of the problem. The full code can be found below the stars-line but I have copied the portion of interest here:

*****************************************************************************

"";" Long-run parameter estimates and s.e.";

bt1 = ththat1[3:k2+1]; phi1=ththat1[2];
p1 = (eye(k2-1)~(bt1./(1-phi1)))./(1-phi1);
bse1=p1*sb1[2:k2+1,2:k2+1]*p1';

bt2 = tht2hat[2:rows(tht2hat)]; phi2=tht2hat[1];
p2 = (eye(k2-1)~(bt2/(1-phi2)))/(1-phi2);
bse2=p2*se2*p2';

bt = bt1|bt2; bse = sqrt(diag(bse1))|sqrt(diag(bse2));
bt'|bse';

*****************************************************************************

Thanks for your help.

Ozey

 

 

new; format /rd /m1 8,4; output file= shoot.out on;
tstart = date;
fname="mydata.xlsx";
range="a2:d821";
sheet=1;
bidata=spreadSheetReadM(fname,range,sheet);
t = 20;
nt = rows(bidata);
n = nt/t;

y = bidata[.,1]; @ real gdp growth rate @
iiq = bidata[.,2]; @ index of institutional quality @
iis = bidata[.,3]; @ index of infrastructure stock @
pcp = bidata[.,4]; @ physical capital @
hcp = bidata[.,5]; @ human capital @
top = bidata[.,6]; @ trade openness @
goc = bidata[.,7]; @ government consumption @

/***************************** Controls *******************************/

qn = 200; @ number of quantiles to examine @
_trim_1 = .15; @ percentage to trim before search, single @
_con = 0; if _con ==0; "no constant"; elseif _con==1;
"constant in the upper regime";endif;
ns = 200; @ iter num in averaging (1st uses analytic formula) @
nb = 0; @ bootstrap iteration @
jmy = 1; @ number of lags of y used in IV ( jm >= 1) @
jmx = 1; @ number of lags of x used in IV ( jm >= 1) @
jn = 1; @ initial lag used in IV is t-jn-2 (jn >=0) @
t0 = jn+3; @ initial time in GMM ( >= jn+3) @
jj=1000; @ iteration for linearity test p-value @

x = iiq~iis~pcp~hcp~top~goc; @ set regressors @
// tt= ones(n,1).*.eye(t); tt= tt[.,t0+1:t]; // time dummies. SET _con =0
// x = x~tt; with dummies we need more IV's as they are already part of IVs

q = iiq;

xx= ones(nt,1)~iis~hcp; @ variables to be used as IV: FIRST COLUMN MUST BE 1's@

"threshold variable = iiq ";
"qn~ns~jm~jn~t0 : " qn~ns~jmy~jn~t0;

/*************************************************************************/

w = setiv(n,t,y,xx,t0,jmy,jmx,jn);" number of moments: " cols(w);

dd = unique(q,1);
qnt1 = qn*_trim_1;
sq = seqa(_trim_1,1/qn,qn-2*qnt1+1);
qq1 = dd[floor(sq*rows(dd))]; @ grid for threshold values @
qn1 = rows(qq1);

ly = vec(lagn(reshape(y,n,t)',1)); // ntx1
dy = y-ly; // ntx1
z2 = ly;
k=cols(x); z2 = ly~x;
z1 = z2;
if _con == 1; z1= (z2~ones(n*t,1)); endif;
k1 = cols(z1); k2 = cols(z2);

dz0 = {} ; @ stack of dz's across threshold values @
for i (1,qn1,1);

z = z2~(z1.*(q.> qq1[i])); // ntx[(1+k)+(1+k)+1]
lz = lagnmatrix(n,t,z,1); // ntx[(1+k)+(1+k)]
dz = z-lz; // ntx[(1+k)+(1+k)]
tmp = packr(dy~dz); // n(t-2)x[1+(1+k)+(1+k)]
tmp = trimrmatrix(n,t-2,tmp,t0-3,0);
dz = tmp[.,2:cols(tmp)]; // n(t-t0+1)x(1+k+1+1+k)
dz0=dz0|dz;

endfor; dy = tmp[.,1];

// ************************** 1st step GMM

zst = rndn(cols(w),jj); // random numbers for linearity test
thtsample = {}; Jsample={};
wnst_A = zeros(qn1,jj); wn_A = zeros(qn1,1);
colp = zeros(ns,1);

for iter (1,ns,1);

if iter == 1; iv = weight1(n,t,w,t0); else;
e = rndn(t-1,n); // (t-1)xn, simulating homogeneous case
de = vec(e[2:t-1,.]-e[1:t-2,.]); // n(t-2)x1
iv = weight(n,t,de,w,t0);
endif;

col1={};
for i (1,qn1,1);

dz = dz0[(i-1)*n*(t-t0+1)+1:i*n*(t-t0+1),.];
wdz = (w'dz)/n; // mx(1+k+1+1+k)
wdy = (w'dy)/n; // mx1
trap 1; iwdz = inv(wdz'*iv*wdz); trap 0;
if scalerr(iwdz); iwdz = invswp(wdz'*iv*wdz);
" error in the 1st step " i; endif;

b = iwdz*(wdz'*iv*wdy); // (1+k+1+k)x1
wde = wdy-wdz*b; // mx1
s = wde'*iv*wde; // scalar
dehat = dy-dz*b; // n(t-2)x1

col1=col1|(s~qq1[i]~b'~dehat');

endfor;

col1hat=sortc(col1,1);
dehat = col1hat[1,2+rows(b)+1:cols(col1)]';

// ************************** 2nd step GMM

iv = weight(n,t,dehat,w,t0);
wnst=zeros(qn1,jj);
wn = zeros(qn1,1);

col2={};
for i (1,qn1,1);

dz = dz0[(i-1)*n*(t-t0+1)+1:i*n*(t-t0+1),.];
wdz = (w'dz)/n; // mx(1+k+1+1+k)
wdy = (w'dy)/n; // mx1
trap 1; iwdz = inv(wdz'*iv*wdz); trap 0;
if scalerr(iwdz); iwdz = invswp(wdz'*iv*wdz);
" error in the 2nd step " i; endif;
b = iwdz*(wdz'*iv*wdy); // (1+k+1+k)x1

wde = wdy-wdz*b; // mx1
s = wde'*iv*wde; // scalar
dehat = dy-dz*b; // n(t-2)x1

col2=col2|(s~qq1[i]~b'~dehat');

//-- supW statistic
p1 = zeros(k1,k2)~eye(k1);
ivi = weight(n,t,col1[i,2+rows(b)+1:cols(col1)]',w,t0);
iv2 = chol(ivi);
trap 1; iwdzi = inv(wdz'*ivi*wdz); trap 0;
if scalerr(iwdzi); iwdzi = invswp(wdz'*ivi*wdz);
" error in the 2nd step " i; endif;

idti = inv(p1*iwdzi*p1');
wn[i] = n*(p1*b)'idti*(p1*b);

for j (1,jj,1);
wnst[i,j] = zst[.,j]'iv2*wdz*iwdzi*p1'idti*p1*iwdzi*wdz'iv2'zst[.,j];
endfor;

endfor;

col2hat=sortc(col2,1);
dehat = col2hat[1,2+rows(b)+1:cols(col2hat)]';

thtsample = thtsample|col2hat[1,2:2+rows(b)];
Jsample = Jsample|(n*col2hat[1,1]);

colp[iter] = 1-counts(maxc(wnst),maxc(wn))/jj; //linearity test
wnst_A = wnst_A + wnst; wn_A = wn_A + wn;

endfor;

ththat = meanc(thtsample); // averaged estimate, (1+k+1+k)x1
ththat1= thtsample[1,.]'; // estimate from analytic Weight

Jhat = meanc(Jsample); // averaged J statistic
Jhat1= Jsample[1]; // J stat from analytic Weight

ltp_A = 1-counts(maxc(wnst_A),maxc(wn_A))/jj;
"";" linearity test from analytic : " colp[1];
" linearity test from averaging: " ltp_A; "";

""; "/------ output based on analytic variance formula -----/";"";;
{dum1,dum2,dum3,dum4}=printout(n,t,t0,w,dy,q,ththat1,z1,z2,Jhat1,_con);
//dum1: s.e.; dum2: upper regime estimates and s.e.;
//dum3: long-run parameter estimates and s.e. dum4: Jstat and p-value

if ns > 1; ""; "/------ in case of averaging -----/";"";;
{dum1,dum2,dum3,dum4}=printout(n,t,t0,w,dy,q,ththat,z1,z2,Jhat,_con);
endif;

tend = date;et = etstr(ethsec(tstart,tend)); "";" elapsed time : " et;
end;

/*****************************************************
input :

t0: initial time for the moment condition
jm's: maximum lags to be used as an IV
x: should contain ones in the first column

*******************************************************/

proc setiv(n,t,y,x,t0,jmy,jmx,jn);
local er,er1,j,k,m,i,zyi,zxi,wi,sr,w,sj,mj;

k = cols(x); // x[.,1]=1
w = {};
for i (1,n,1);
zyi = y[(i-1)*(t)+1:i*(t)]; // (t)x1
zxi = x[(i-1)*(t)+1:i*(t),.]; // (t)xk
wi = {}; // (t-t0+1)xm
er = 0; er1=0; // scalar

for j (1,t-t0+1,1);
sj = j+t0-jn-3;
mj = 1+ minc(jmy|sj);
wi = wi~zeros(t-t0+1,mj);
sr = er+1; // scalar
er1 = er+ mj; // scalar
er = er1;
wi[j,sr:er1]=1~zyi[maxc(1|sj-jmy+1):sj]';

if k > 1;
mj = (-1+k)*minc(jmx|sj);
wi = wi~zeros(t-t0+1,mj);
er = er1+ mj; // scalar
wi[j,er1+1:er]= vec(zxi[maxc(1|sj-jmx+1):sj,2:k])';
endif;
endfor;

w = w|wi; // (t-t0+1)xm
endfor;
retp(w); endp;

/* input:
de: residuals
w : IV's
*/

proc weight(n,t,de,w,t0); // with given residuals
local v,vbar,i,dei,wi,iv;

v = 0; // scalar
vbar = 0;

for i (1,n,1);
dei = de[(i-1)*(t-t0+1)+1:i*(t-t0+1),.]; // (t-2)x1
wi = w[(i-1)*(t-t0+1)+1:i*(t-t0+1),.]'dei; // m x 1
v = v + wi*wi'; // mxm
vbar = vbar + wi; // mx1

endfor;
if rank(v/n-(vbar/n)*(vbar/n)') == cols(w);
iv = inv(v/n-(vbar/n)*(vbar/n)'); // mxm
else; "error 2";
iv = invswp(v/n-(vbar/n)*(vbar/n)');
endif;

retp(iv); endp;

proc weight1(n,t,w,t0); //analytic
local m,v,i,j,wi,wj,iv;

m = cols(w); // number of moments
v = zeros(m+2,m+2);

for i (1,n,1);
wi = zeros(t-t0+3,m+2);
wi[2:t-t0+2,2:m+1] = w[(i-1)*(t-t0+1)+1:i*(t-t0+1),.];

for j (1,t-t0+2,1);
wj = wi[j,.]'-wi[j+1,.]';
v = v + wj*wj'; // mxm
endfor;
endfor;

v = v[2:(m+1),2:(m+1)];
if rank(v/n) == cols(w);
iv = inv(v/n); // mxm
else; "error";
iv = invswp(v/n);
endif;

retp(iv); endp;

// ===========================================================================
proc (1) = lagnmatrix(n,t,x,p);
local k,lx,i;
k = cols(x);
lx = zeros(n*t,k);
for i (1,k,1);
lx[.,i] = vec(lagn(reshape(x[.,i],n,t)',p)); // ntx1
endfor;

retp(lx); endp;

proc (1) = trimrmatrix(n,t,x,p,q);
local k,y,i;
k = cols(x);
y = zeros(n*(t-p-q),k);
for i (1,k,1);
y[.,i] = vec(trimr(reshape(x[.,i],n,t)',p,q));
endfor;
retp(y); endp;

// ===========================================================================
//outputs:
//dum1: s.e.; dum2: upper regime estimates and s.e.;
//dum3: long-run parameter estimates and s.e. dum4: Jstat and p-value

proc (4) = printout(n,t,t0,w,dy,q,ththat1,z1,z2,Jhat1,_con);
local dumm1,z,lz,dz,tmp,dehat1,iv,Gb,h,kn,s,ds,i,si,Gr,
sb1,sg1,se1,k2,p,tht2hat,se2,bt1,phi1,p1,bse1,bt2,phi2,p2,bse2,bt,bse;
// *********************************** COVARIANCE MATRIX
dumm1 = (q.> ththat1[1]); " proportion of upper regime: " meanc(dumm1);"";
z = z2~(z1.*dumm1); // ntx[(1+k)+(1+k)+1]
lz = lagnmatrix(n,t,z,1); // ntx[(1+k)+(1+k)]
dz = z-lz; // ntx[(1+k)+(1+k)]
tmp = packr(dz); // n(t-2)x[1+(1+k)+(1+k)]
tmp = trimrmatrix(n,t-2,tmp,t0-3,0);
dz = tmp[.,1:cols(tmp)]; // n(t-t0+1)x(1+k+1+1+k)
dehat1 = dy-dz*ththat1[2:rows(ththat1)]; // n(t-t0+1)x1

iv = weight(n,t,dehat1,w,t0);

Gb = -w'dz/n; // mx[(1+k)+(1+k)]

// *********************************** ds: n(t-2)x(1+k)
h = 1.06*stdc(packr(q))*n^(-0.2); // scalar, bandwidth
kn = exp(-(((ththat1[1]-q)/h).^2)/2)/sqrt(2*pi); // ntx1
s = z1.*kn; // ntx(1+k)
ds = zeros(n*(t-t0+1),cols(s)); // n(t-2)x(1+k)
for i (1,cols(s),1);
si = reshape(s[.,i],n,t)'; // txn
ds[.,i] = vec(trimr(si,t0-1,0)-trimr(si,t0-2,1)); // (t-2)x(1+k)
endfor;

Gr = -(w'ds/(n*h))*ththat1[3+k:rows(ththat1)]; // mx1

// *********************************** STANDARD ERROR
sb1 = inv(Gb'*iv*Gb-(Gb'*iv*Gr)*inv(Gr'*iv*Gr)*(Gr'*iv*Gb))/n;
sg1 = sqrt(1/(n*(Gr'*iv*Gr-(Gr'*iv*Gb)*inv(Gb'*iv*Gb)*(Gb'*iv*Gr))));
se1 = sg1|sqrt(diag(sb1));

" estimates and s.e. (threshold~lower regime (lag y,x)~delta (1,lag y,x)";
ththat1'|se1';

"";" upper regime estimates and s.e. (lag y,x)";

k2 = cols(z2);
if _con ==1; p = eye(k2)~zeros(k2,1)~eye(k2);
else; p = eye(k2)~eye(k2); endif;
tht2hat = (zeros(k2,1)~p)*ththat1;
se2 = p*sb1*p';
tht2hat'|sqrt(diag(se2))';

"";" Long-run parameter estimates and s.e.";

bt1 = ththat1[3:k2+1]; phi1=ththat1[2];
p1 = (eye(k2-1)~(bt1./(1-phi1)))./(1-phi1);
bse1=p1*sb1[2:k2+1,2:k2+1]*p1';

bt2 = tht2hat[2:rows(tht2hat)]; phi2=tht2hat[1];
p2 = (eye(k2-1)~(bt2/(1-phi2)))/(1-phi2);
bse2=p2*se2*p2';

bt = bt1|bt2; bse = sqrt(diag(bse1))|sqrt(diag(bse2));
bt'|bse';

"";" Overidentification J-statistic (p-value) " ;
Jhat1~cdfchic(Jhat1,cols(W)-rows(ththat1));

retp(se1,tht2hat~sqrt(diag(se2)),bt~bse,Jhat1~cdfchic(Jhat1,cols(W)-rows(ththat1)));
endp;

 

 

Your Answer


You must login to post answers.

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.