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;