Dear All,
in my research I need to replicate the first stage of the spivregress command. Unfortunately the routines do not offer this option. Following user manual, I came out with the following code in Mata:
It seems to run smoothly (I hope it is correct: I am new in Mata) until the very end. When I reach the lines in red I get the following error message:
I am running out of ideas. Do you have any suggestion? I would be very grateful for any suggestion you may want to provide.
Best regards,
Dario
in my research I need to replicate the first stage of the spivregress command. Unfortunately the routines do not offer this option. Following user manual, I came out with the following code in Mata:
Code:
********************************************************************************
* Replicate spivregress first-stage (e(delta_2sls)) and final GS2SLS (e(b))
* Manual first stage (2SLS) and manual second stage (GS2SLS) for SARAR(1,1)
********************************************************************************
clear all
version 18.0
set matsize 11000
mata:
mata clear
mata set matastrict on
end
use dataset.dta, clear
replace literacy1880 = 46.77 if NUTS_ID=="UKN" /* Add the literacy rate for Norther Ireland */
replace literacy1880 = 24.33 if NUTS_ID=="ES53" /* Add the literacy rate for Baleares */
spset, modify coordsys(latlong, kilometers)
spmatrix create idistance dW1500, vtruncate(1/1500) replace
spmatrix summarize dW1500
********************************************************************************
* 0) USER SETTINGS
********************************************************************************
local yvar yp9500
local endog pc_culture // <-- confirm spelling matches dataset
local ziv pc_institutions literacy1880
local Xlist urb_rate1850 school country1 country2 country3 country4 country5 country7 country8 country15
local Wname dW1500 // spatial matrix name
local IMP 2 // impower used in spivregress
local TOL 1e-10 // tolerance for collinearity filtering
********************************************************************************
* 1) Spatial setup & spivregress
********************************************************************************
spset
spmatrix dir
* Bring W into Mata
spmatrix matafromsp W id = `Wname'
* Run spivregress (heteroskedastic SARAR(1,1), no constant)
spivregress `yvar' (`endog' = `ziv') `Xlist', ///
heteroskedastic dvarlag(`Wname') errorlag(`Wname') ///
noconstant impower(`IMP')
* Store results
matrix SPIV_delta2sls = e(delta_2sls)
matrix SPIV_b = e(b)
scalar SPIV_rho = _b[`Wname':e.`yvar'] // SAR error parameter
* Detect Wy coefficient name for comparison
local ebnames : colfullnames e(b)
local wyname
foreach nm of local ebnames {
if (strpos("`nm'", "`Wname'") & strpos("`nm'", "`yvar'")) local wyname `nm'
}
if "`wyname'" == "" {
foreach nm of local ebnames {
if (strpos("`nm'", " W ") & strpos("`nm'", "`yvar'")) local wyname `nm'
}
}
if "`wyname'" != "" {
scalar BLag = _b["`Wname':`yvar'"]
di as res "Detected Wy coefficient name: `wyname' | value: " %9.6f BLag
}
di as res "rho_hat (SAR error): " %9.6f SPIV_rho
********************************************************************************
* 2-3) MATA — Definitions and Execution
********************************************************************************
mata:
struct colfilter {
real matrix X
real colvector keep
}
struct firststage {
real matrix Z
real matrix H1
real colvector theta2sls
real scalar rankH1
real scalar rankZ
}
struct gs2sls {
real colvector beta_gs2sls
real scalar rankH2
}
struct colfilter drop_dep_cols(real matrix X, real scalar tol) {
struct colfilter CF;
CF.keep = J(cols(X), 1, 0);
real matrix kept;
kept = J(rows(X), 0, .);
real matrix basis;
basis = J(0, 0, .);
real scalar cur_rank;
cur_rank = 0;
real scalar j;
real matrix trial;
real matrix U, V;
real colvector s;
real scalar r2;
for (j = 1; j <= cols(X); j++) {
trial = (basis, X[, j]);
svd(trial, U, s, V);
r2 = sum(s :> tol);
if (r2 > cur_rank) {
kept = kept , X[, j];
CF.keep[j] = 1;
basis = trial;
cur_rank = r2;
}
}
CF.X = kept;
return(CF);
}
struct firststage spiv_first_stage_singleW(real colvector y, real matrix X, real matrix Yend, real matrix Ziv, real matrix W, real scalar impower, real scalar tol) {
real matrix X53;
real matrix H;
real matrix curW;
real scalar r;
real matrix H1;
real colvector Wy;
real matrix HtH;
real matrix invHtH;
real matrix PH1;
real matrix Ztil;
real matrix A;
real colvector b;
if (impower <= 0) impower = 2;
if (tol <= 0) tol = 1e-10;
X53 = X , Ziv;
H = X53;
curW = W;
for (r = 1; r <= impower; r++) {
H = H , (curW * X53);
curW = curW * W;
}
struct colfilter CF2;
CF2 = drop_dep_cols(H, tol);
H1 = CF2.X;
Wy = W * y;
real matrix Z;
Z = Yend , X , Wy;
HtH = quadcross(H1, H1);
invHtH = invsym(HtH);
PH1 = H1 * invHtH * H1';
Ztil = PH1 * Z;
A = quadcross(Ztil, Z);
b = Ztil' * y;
struct firststage FS;
FS.theta2sls = invsym(A) * b;
FS.Z = Z;
FS.H1 = H1;
FS.rankH1 = rank(H1);
FS.rankZ = rank(Z);
return(FS);
}
struct gs2sls spiv_second_stage_gs2sls_singleW(real colvector y, real matrix Z, real matrix H1, real matrix W, real scalar rho_hat, real scalar tol) {
if (tol <= 0) tol = 1e-10;
real matrix H2full;
H2full = H1 , (W * H1);
struct colfilter CF3;
CF3 = drop_dep_cols(H2full, tol);
real matrix H2;
H2 = CF3.X;
real colvector Wy;
Wy = W * y;
real colvector ystar;
ystar = y - rho_hat * Wy;
real matrix Zstar;
Zstar = Z - rho_hat * (W * Z);
real matrix HtH;
HtH = quadcross(H2, H2);
real matrix invHtH;
invHtH = invsym(HtH);
real matrix PH2;
PH2 = H2 * invHtH * H2';
real matrix Ztil;
Ztil = PH2 * Zstar;
real matrix A1;
A1 = quadcross(Ztil, Zstar);
real matrix b1;
b1 = Ztil' * ystar;
struct gs2sls OUT;
OUT.beta_gs2sls = invsym(A1) * b1;
OUT.rankH2 = rank(H2);
return(OUT);
}
struct manual{
real colvector y1
real matrix X1
real matrix Yend1
real matrix Ziv1
real matrix W
real scalar imp
real scalar tlr
}
y1 = st_data(., tokens("`yvar'"));
X1 = st_data(., tokens("`Xlist'"));
Yend1 = st_data(., tokens("`endog'"));
Ziv1 = st_data(., tokens("`ziv'"));
imp = strtoreal("`IMP'");
tlr = strtoreal("`TOL'");
FS = spiv_first_stage_singleW(y1, X1, Yend1, Ziv1, W, imp, tlr);
st_matrix("MAN_delta2sls", FS.theta2sls);
real scalar rho = st_numscalar("SPIV_rho");
struct gs2sls G = spiv_second_stage_gs2sls_singleW(y1, FS.Z, FS.H1, W, rho, tlr);
st_matrix("MAN_beta_gs2sls", G.beta_gs2sls);
end
********************************************************************************
* 4) Compare results
********************************************************************************
di as txt "==> Comparing first-stage (Stata vs Manual):"
mat list SPIV_delta2sls, format(%9.6f)
mat list MAN_delta2sls , format(%9.6f)
di as txt "==> Comparing final-stage (GS2SLS) coefficients:"
tempname Bsubset
matrix `Bsubset' = J(colsof(MAN_beta_gs2sls), 1, .)
local pos 1
matrix `Bsubset'[`pos',1] = _b[`endog']
local ++pos
foreach v of local Xlist {
matrix `Bsubset'[`pos',1] = _b[`v']
local ++pos
}
if "`wyname'" != "" {
matrix `Bsubset'[`pos',1] = BLag
}
mat list `Bsubset', format(%9.6f) title("Stata e(b) subset [endog, X, W*y]")
mat list MAN_beta_gs2sls, format(%9.6f) title("Manual GS2SLS [endog, X, W*y]")
Code:
: FS = spiv_first_stage_singleW(y1, X1, Yend1, Ziv1, W, imp, tlr);
drop_dep_cols(): 3259 nonstruct found where struct required
spiv_first_stage_singleW(): - function returned error
<istmt>: - function returned error
(6 lines skipped)
--------------------------------------------------------------------------------------------
r(3259);
end of do-file
r(3259);
Best regards,
Dario
Comment