# delimit; clear; clear matrix; clear mata; set mem 100m; set matsize 800; cd "c:\Research\Dynamic\"; capture log close; log using dynamic02_PSID.log, replace; scalar T=12; // Does not count the initial condition; use sw-data.dta; qui sum id if year==1; local nobs=r(N); scalar N=`nobs'; ************ DO FD-IV ESTIMATION **************; replace lnw=. if s==0; tsset id year; gen lnw_lag=l.lnw; replace lnw_lag=lnw80 if year==1; gen lnw_lag2=l.lnw_lag; gen dlnw=lnw-lnw_lag; gen dlnw_lag=d.lnw_lag; gen dagesq=d.agesq; gen d=(lnw~=. & lnw_lag~=. & lnw_lag2~=.); local t=2; while `t'<=T {; qui probit d age agesq educ children children_lag1 children_lag2 if year==`t'; qui predict xb, xb; qui gen lam`t'= 0; qui replace lam`t'= normalden(xb)/normal(xb)*d if year==`t'; qui drop xb; local t=`t'+1; }; tab year, gen(td); ivreg dlnw dagesq td* (dlnw_lag=lnw_lag2), robust cluster(id); scalar rho_FDIV_nlam=_b[dlnw_lag]; scalar b_agesq_FDIV_nlam=_b[dagesq]; scalar se_rho_FDIV_nlam=_se[dlnw_lag]; scalar se_agesq_FDIV_nlam=_se[dagesq]; ivreg dlnw dagesq lam* td* (dlnw_lag=lnw_lag2), robust cluster(id); test lam2 lam3 lam4 lam5 lam6 lam7 lam8 lam9 lam10 lam11 lam12; scalar rho_FDIV_lam=_b[dlnw_lag]; scalar b_agesq_FDIV_lam=_b[dagesq]; scalar se_rho_FDIV_lam=_se[dlnw_lag]; scalar se_agesq_FDIV_lam=_se[dagesq]; mat b=e(b); mat V=e(V)'; mat b1=b[1,3..13]; mat V1=V[3..13,3..13]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2_FDIV=chi2[1,1]; scalar pv_ch2_FDIV=chi2tail(11,ch2_FDIV); scalar di ch2_FDIV; scalar di pv_ch2_FDIV; more; drop d lam* td* children_lag1 children_lag2; ***************************************************************************; capture program drop nldyn; // NLS on main equation, no lambda; program nldyn; version 11; syntax varlist if, at(name); local mysslhs: word 1 of `varlist'; local myinicd: word 2 of `varlist'; local myssrhs; local nv=3; while `nv'<=NV {; local myssrhs`nv': word `nv' of `varlist'; local myssrhs `myssrhs' `myssrhs`nv''; local nv=`nv'+1; }; // Obtain parameters on 2nd step time-constant (x1,...xt,z1,...,zt) variables // and compute the total contribution of these variables to be // multiplied by rho function later, when obtaining moment conditions; tempvar ssalltc; qui gen double `ssalltc'=0 `if'; local k=0; foreach var of varlist `myssrhs' {; local k=`k'+1; qui replace `ssalltc'= `ssalltc' + `var'*`at'[1,`k'] `if'; }; scalar nsstcrhs=`k'; // Excludes the intercept; * noisily scalar di nsstcrhs; * noisily di "2nd step cum eff. of tc vars, ssalltc"; * noisily sum `ssalltc'; // Obtain parameters on second-step time-varying (x) variables // (following time-constant second-step parameters); foreach var of varlist $sstvvars {; local k=`k'+1; scalar beta`var'=`at'[1,`k']; * noisily scalar di beta`var'; }; scalar nssrhs=`k'; // Excludes a constant and initial condition; * noisily scalar di nssrhs; // Obtain parameter eta11...eta1T (following 2nd step x & z variables); local t=1; while `t'<=T {; scalar eta1`t'=`at'[1,nssrhs+`t']; local t=`t'+1; }; // Obtain parameter gamma (following 2nd step x & z variables and eta1); scalar gamma=`at'[1,nssrhs+T+1]; // Obtain parameter rho (following 2nd step x & z variables, eta1 and gamma); scalar rho=`at'[1,nssrhs+T+2]; // Generate starting values for time-varying variables // to be used later for computing moment conditions at the 2nd step; foreach var of varlist $sstvvars {; // Accum. second-step time-varying X; tempvar cum`var'0; qui gen double `cum`var'0'=0 `if'; }; tempvar cumcons0; qui gen double `cumcons0'=0 `if'; // Contribution of all time-varying variables; * tempvar ssalltv; capture drop ssalltv; qui gen double ssalltv=0 `if'; tempvar lam_all; qui gen double `lam_all'=0 `if'; //Start the cycle over time - compute 2nd step conditional mean depending on t; local t=1; while `t'<=T {; * noisily di "Time period t = " `t'; // Add up all cumulative contributions of second-step time-varying X // at time t; local t_lag=`t'-1; foreach var of varlist $sstvvars {; tempvar cum`var'`t'; qui gen double `cum`var'`t''=rho*`cum`var'`t_lag''+`var'`t' `if'; qui replace ssalltv=ssalltv+`cum`var'`t''*beta`var' if year==`t'; }; // Add the cumulative effect of time-specific intercept tempvar cumcons`t'; qui gen double `cumcons`t''=rho*`cumcons`t_lag''+eta1`t' `if'; qui replace ssalltv=ssalltv+`cumcons`t'' if year==`t'; local t=`t'+1; }; // Generate 2-step moment condition for time period t; qui replace `mysslhs'= s*( (rho^year)*`myinicd' + ssalltv + ((1-rho^year)/(1-rho))*(`ssalltc' + gamma*`myinicd')) `if'; end; ***************************************************************************; ***************************************************************************; capture program drop nldselect; // NLS on main equation, estimated lambda; program nldselect; version 11; syntax varlist if, at(name); local mysslhs: word 1 of `varlist'; local myinicd: word 2 of `varlist'; local myssrhs; local nv=3; while `nv'<=NV {; local myssrhs`nv': word `nv' of `varlist'; local myssrhs `myssrhs' `myssrhs`nv''; local nv=`nv'+1; }; // Obtain parameters on 2nd step time-constant (x1,...xt,z1,...,zt) variables // and compute the total contribution of these variables to be // multiplied by rho function later, when obtaining moment conditions; tempvar ssalltc; qui gen double `ssalltc'=0 `if'; local k=0; foreach var of varlist `myssrhs' {; local k=`k'+1; qui replace `ssalltc'= `ssalltc' + `var'*`at'[1,`k'] `if'; }; scalar nsstcrhs=`k'; // Excludes the intercept; * noisily scalar di nsstcrhs; * noisily di "2nd step cum eff. of tc vars, ssalltc"; * noisily sum `ssalltc'; // Obtain parameters on second-step time-varying (x) variables // (following time-constant second-step parameters); foreach var of varlist $sstvvars {; local k=`k'+1; scalar beta`var'=`at'[1,`k']; * noisily scalar di beta`var'; }; scalar nssrhs=`k'; // Excludes a constant and initial condition; * noisily scalar di nssrhs; // Obtain parameter eta11...eta1T (following 2nd step x & z variables); local t=1; while `t'<=T {; scalar eta1`t'=`at'[1,nssrhs+`t']; local t=`t'+1; }; // Obtain parameter gamma (following 2nd step x & z variables and eta1); scalar gamma=`at'[1,nssrhs+T+1]; // Obtain parameter rho (following 2nd step x & z variables, eta1 and gamma); scalar rho=`at'[1,nssrhs+T+2]; // Generate starting values for time-varying variables // to be used later for computing moment conditions at the 2nd step; foreach var of varlist $sstvvars {; // Accum. second-step time-varying X; tempvar cum`var'0; qui gen double `cum`var'0'=0 `if'; }; tempvar cumcons0; qui gen double `cumcons0'=0 `if'; // Contribution of all time-varying variables; * tempvar ssalltv; capture drop ssalltv; qui gen double ssalltv=0 `if'; tempvar lam_all; qui gen double `lam_all'=0 `if'; //Start the cycle over time - compute 2nd step conditional mean depending on t; local t=1; while `t'<=T {; * noisily di "Time period t = " `t'; // Add up all cumulative contributions of second-step time-varying X // at time t; local t_lag=`t'-1; foreach var of varlist $sstvvars {; tempvar cum`var'`t'; qui gen double `cum`var'`t''=rho*`cum`var'`t_lag''+`var'`t' `if'; qui replace ssalltv=ssalltv+`cum`var'`t''*beta`var' if year==`t'; }; // Add the cumulative effect of time-specific intercept tempvar cumcons`t'; qui gen double `cumcons`t''=rho*`cumcons`t_lag''+eta1`t' `if'; qui replace ssalltv=ssalltv+`cumcons`t'' if year==`t'; // Add up contribution of all lambdas; qui replace `lam_all'=`lam_all'+lam`t'*`at'[1,nssrhs+T+2+`t']; local t=`t'+1; }; // Generate 2-step moment condition for time period t; qui replace `mysslhs'= s*( (rho^year)*`myinicd' + ssalltv + ((1-rho^year)/(1-rho))*(`ssalltc' + gamma*`myinicd') + `lam_all') `if'; qui replace g = `mysslhs' `if'; * noisily sum mm; * noisily scalar di betax; * noisily scalar di rho; * noisily scalar di gamma; end; ***************************************************************************; qui gen g=.; sort id year; local t=1; while `t'<=T {; foreach var of varlist age agesq children {; by id: gen `var'`t'=`var'[`t']; }; local t=`t'+1; }; rename age tvage; rename agesq tvagesq; rename children tvchildren; qui gen age=.; qui gen agesq=.; global sstvvars age agesq; local t=1; while `t'<=T {; qui probit s tvage tvagesq educ lnw80 children* if year==`t'; qui predict xb, xb; qui gen lam`t'= 0; qui replace lam`t'= normalden(xb)/normal(xb) if year==`t' & s==1; qui drop xb; local t=`t'+1; }; * DO NLS W/O SELECTION CORRECTION; mat ivals=J(1,29,0); scalar NV=15; // NUMBER OF VARIABLES IN UNOBSERVED EFFECT + DEPVAR + INICOND; nl dyn @ lnw lnw80 educ children*, nparameters(29) initial(ivals) vce(cluster id); scalar b_educ_NLS_nlam=_b[/b1]; scalar b_age_NLS_nlam=_b[/b14]; scalar b_agesq_NLS_nlam=_b[/b15]; scalar rho_NLS_nlam=_b[/b29]; scalar se_educ_NLS_nlam=_se[/b1]; scalar se_age_NLS_nlam=_se[/b14]; scalar se_agesq_NLS_nlam=_se[/b15]; scalar se_rho_NLS_nlam=_se[/b29]; * DO NLS WITH SELECTION CORRECTION; mat ivals=J(1,41,0); scalar NV=15; // NUMBER OF VARIABLES IN UNOBSERVED EFFECT + DEPVAR + INICOND; nl dselect @ lnw lnw80 educ children*, nparameters(41) initial(ivals) vce(cluster id) leave; predict ehat, res; scalar b_educ_NLS_lam=_b[/b1]; scalar b_age_NLS_lam=_b[/b14]; scalar b_agesq_NLS_lam=_b[/b15]; scalar rho_NLS_lam=_b[/b29]; mat b=e(b); mat VR=e(V); mat b1=b[1,30..41]; mat V1=VR[30..41,30..41]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2_NLS=chi2[1,1]; scalar pv_ch2_NLS=chi2tail(12,ch2_NLS); scalar di ch2_NLS; scalar di pv_ch2_NLS; gen res=(lnw-g); gen res2=(lnw-g)^2; sum ehat res res2; *OBTAIN COEFFICIENTS FOR LAMBDA TERMS FROM THE REGRESSION; local i=30; local j=1; while `i'<=41 {; scalar psi`j'=_b[/b`i']; scalar di psi`j'; local i=`i'+1; local j=`j'+1; }; *************************************************************; * OBTAINING CLUSTER-ROBUST STANDARD ERRORS MANUALLY; foreach var of varlist b1-b41 {; qui replace `var'=0 if lnw==.|lnw==0; }; local varli2; /* Here accumulate sum (over t) of derivatives of E(y) wrt parameter x residual */ foreach var of varlist b1-b41 {; qui gen eh`var'=`var'*ehat; qui egen t`var' = sum(eh`var'), by(id); local varli2 `varli2' t`var'; }; sum `varli2'; ********************************************************************; * SCALAR g BELOW IS THE NUMBER OF INDIVIDUALS IN THE SELECTED SAMPLE ********************************************************************; sort id year; by id: gen num=_n; gen countid=(num==1); sum countid; scalar g=r(sum); drop countid; mat accum MM=b1-b41, nocons; mat accum MEEM=`varli2' if num==1, nocons; mat VmyR=syminv(MM)*MEEM*syminv(MM); *DF ADJUSTMENT, NOT USED: (e(N)-1)*g/((g-1)*(e(N)-e(df_m))); scalar se_educ=sqrt(VmyR[1,1]); scalar se_age=sqrt(VmyR[14,14]); scalar se_agesq=sqrt(VmyR[15,15]); scalar se_rho=sqrt(VmyR[29,29]); scalar di se_educ; scalar di se_age; scalar di se_agesq; scalar di se_rho; more; ******************************************************************** * OBTAIN STANDARD ERRORS CORRECTED FOR THE 1-st STEP ESTIMATION ********************************************************************; ******************************************************************** * * GENERATE NEW VARS AND A DIAGONAL MATRIX WITH HESSIANS * * VARS1= GAMMA_t*d(LAMBDA)/d(xb)*VAR BY t, denote them "g" * * VARS2= d(LnL)/d(xb)*VAR BY t - SCORES FOR EACH i, denote them "q" * * DIAG(H)=diag(H_t), t=1,...,T * ********************************************************************; gen cons=1; *****************************************************************; * varli1 BELOW DEFINES THE LIST OF THE VARIABLES USED AS * REGRESSORS AT THE 1st STAGE (PROBIT) * varli3 AND varli4 ARE EMPTY FOR NOW AND WILL BE USED LATER *****************************************************************; local varli1 educ age agesq lnw80 children* cons; local varli3; local varli4; *****************************************************************; * MATRIX H BELOW SHOULD BE A SQUARE MATRIX OF DIMENSION * # REGRESSORS AT THE 1st STAGE (PROBIT) * # TIME PERIODS *****************************************************************; mat H=I(17*12); /* 17 = # vars in probit, 12=T */ local j=0; local t=1; while `t'<=T {; di "Year="`t'; qui probit s educ tvage tvagesq lnw80 children* if year==`t'; predict xb, xb; mat H`t'=e(V); qui gen tempvar1=normalden(xb)/normal(xb) if s==1&year==`t'; qui replace tempvar1=-normalden(xb)/(1-normal(xb)) if s==0&year==`t'; assert lam`t'==tempvar1 if s==1&year==`t'; mat H[17*`j'+1,17*`j'+1]=H`t'; foreach var of varlist educ tvage tvagesq lnw80 children* cons {; qui gen g`var'_`t' = 0 if s==1; qui replace g`var'_`t'= -tempvar1*(tempvar1+xb)*psi`t'*`var' if s==1&year==`t'; qui gen q`var'_`t' = .; qui replace q`var'_`t'= tempvar1*`var' if year==`t'; sort id year; qui by id: replace q`var'_`t'=q`var'_`t'[`t'] if year~=`t'; local varli3 `varli3' g`var'_`t'; local varli4 `varli4' q`var'_`t'; }; local j=`j'+1; drop xb tempvar*; mat drop H`t'; local t=`t'+1; }; ***********************************************************************; * NUMBER OF VARIABLES IN varli2 (WHICH IS THE LIST OF INTERACTION TERMS * (sum over t of each of B1-B41)*) * SHOULD BE EQUAL TO * THE NUMBER OF PARAMERTERS/SCORES AT THE SECOND STAGE, i.e. 41 * * NUMBER OF VARIABLES IN varli4 SHOULD BE EQUAL TO * <# FIRST-STAGE REGRESSORS> * <# TIME PERIODS> ***********************************************************************; mat accum TEMP=`varli2' `varli4' if num==1, nocons; *NOTE: #instruments*ehat (#vars in `varli2')=41; *Number of vars in Q (#vars in `varli4')= 17*12 = 204; ***********************************************************************; * EXTRACT THE UPPER RIGHT CORNER OF THE TEMP MATRIX ***********************************************************************; mat ZEQ=TEMP[1..41,41+1...]; di "Number of rows in ZEQ="rowsof(ZEQ); di "Number of columns in ZEQ="colsof(ZEQ); mat drop TEMP; ***********************************************************************; * NUMBER OF VARIABLES IN varli3 SHOULD BE EQUAL TO * <# FIRST-STAGE REGRESSORS> * <# TIME PERIODS> ***********************************************************************; mat accum TEMP=b1-b41 `varli3', nocons; *NOTE: # vars at 2nd stage=41; *Number of vars in G (#vars in `varli3'= # vars at 1st stage)=17; ***********************************************************************; * EXTRACT THE UPPER RIGHT CORNER OF THE TEMP MATRIX ***********************************************************************; mat ZG=TEMP[1..41,41+1...]; di "Number of rows in ZG="rowsof(ZG); di "Number of columns in ZG="colsof(ZG); mat drop TEMP; *********************************************************************; * IF EVERYTHING ABOVE WAS DONE CORRECTLY, THE REMAINING COMPUTATIONS * BELOW SHOULD FOLLOW AUTOMATICALLY *********************************************************************; mat TERM2=ZEQ*H*ZG'; mat accum QQ=`varli4' if num==1, nocons; mat TERM4=ZG*H*QQ*H*ZG'; mat C=(MEEM-TERM2-TERM2'+TERM4); mat VCR=inv(MM)*C*inv(MM); *RESULTS WITH ROBUST V-C MATRIX THAT ALSO TAKES INTO ACCOUNT THE FIRST-STEP ESTIMATION; scalar se_educ_NLS_lam=sqrt(VCR[1,1]); scalar se_age_NLS_lam=sqrt(VCR[14,14]); scalar se_agesq_NLS_lam=sqrt(VCR[15,15]); scalar se_rho_NLS_lam=sqrt(VCR[29,29]); scalar di se_educ_NLS_lam; scalar di se_age_NLS_lam; scalar di se_agesq_NLS_lam; scalar di se_rho_NLS_lam; mat b1=b[1,30..41]; mat V1=VCR[30..41,30..41]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2_NLS_lam=chi2[1,1]; scalar pv_ch2_NLS_lam=chi2tail(12,ch2_NLS_lam); scalar di ch2_NLS_lam; scalar di pv_ch2_NLS_lam; drop age* agesq* children*; rename tvage age; rename tvagesq agesq; rename tvchildren children; ********* RESHAPE DATA *****************; replace lnw=0 if s==0; keep id year s lnw lnw80 age agesq educ children; reshape wide s lnw age agesq educ lnw80 children, i(id) j(year); rename educ1 yreduc; rename lnw801 ini_lnw; drop lnw80* educ*; rename yreduc educ; local t=1; while `t'<=T {; qui probit s`t' age`t' agesq`t' educ ini_lnw children*; qui predict xb, xb; qui gen lam`t'= normalden(xb)/normal(xb)*s`t'; qui gen double mm`t'=.; qui drop xb; local t=`t'+1; }; ***************************************************************************; capture program drop gmm_dselect; // Main equation, estimated lambda; program gmm_dselect; version 11; syntax varlist if, at(name) mysslhs(varlist) myinicd(varlist) myssrhs(varlist); local t=1; while `t'<=T {; local sslhs`t': word `t' of `mysslhs'; // second-step LHS; local t=`t'+1; }; // Obtain parameters on 2nd step time-constant (x1,...xt,z1,...,zt - educ) variables // and compute the total contribution of these variables to be // multiplied by rho function later, when obtaining moment conditions; tempvar ssalltc; qui gen double `ssalltc'=0 `if'; local k=0; foreach var of varlist `myssrhs' {; local k=`k'+1; qui replace `ssalltc'= `ssalltc' + `var'*`at'[1,`k'] `if'; }; scalar nsstcrhs=`k'; // Excludes the intercept; * noisily scalar di nsstcrhs; * noisily di "2nd step cum eff. of tc vars, ssalltc"; * noisily sum `ssalltc'; // Obtain parameters on second-step time-varying (x) variables - age, agesq // (following time-constant second-step parameters); foreach var of varlist $sstvvars {; local k=`k'+1; scalar beta`var'=`at'[1,`k']; * noisily scalar di beta`var'; }; scalar nssrhs=`k'; // Excludes a constant and initial condition; * noisily scalar di nssrhs; // Obtain parameter eta11...eta1T (following 2nd step x & z variables); local t=1; while `t'<=T {; scalar eta1`t'=`at'[1,nssrhs+`t']; local t=`t'+1; }; // Obtain parameter gamma (following 2nd step x & z variables and eta11...eta1T); scalar gamma=`at'[1,nssrhs+T+1]; // Obtain parameter rho (following 2nd step x & z variables, eta11...eta1T and gamma); scalar rho=`at'[1,nssrhs+T+2]; // Generate starting values for time-varying variables // to be used later for computing moment conditions at the 2nd step; foreach var of varlist $sstvvars {; // Accum. second-step time-varying X & intercept; tempvar cum`var'0; qui gen double `cum`var'0'=0 `if'; * noisily di "cumulative effect of tv variable" `var' ", t=0"; * noisily sum `cum`var'0'; }; tempvar cumcons0; qui gen double `cumcons0'=0 `if'; //Start the cycle over time - compute 2nd step moments for each t; local t=1; while `t'<=T {; * noisily di "Time period t = " `t'; // Add up all cumulative contributions of second-step time-varying X // at time t; local t_lag=`t'-1; tempvar ssalltv`t'; qui gen double `ssalltv`t''=0 `if'; foreach var of varlist $sstvvars {; tempvar cum`var'`t'; qui gen double `cum`var'`t''=rho*`cum`var'`t_lag''+`var'`t' `if'; qui replace `ssalltv`t''=`ssalltv`t''+`cum`var'`t''*beta`var' `if'; }; // Add the cumulative effect of time-specific intercept tempvar cumcons`t'; qui gen double `cumcons`t''=rho*`cumcons`t_lag''+eta1`t' `if'; qui replace `ssalltv`t''=`ssalltv`t''+`cumcons`t'' `if'; * noisily di "cum. contribution of all 2nd step tv X at time " `t'; * noisily sum `ssalltv`t''; // Generate 2-step moment condition for time period t; local m`t': word `t' of `varlist'; qui replace `m`t''= s`t'*(`sslhs`t'' - (rho^`t')*`myinicd' - `ssalltv`t'' - ((1-rho^`t')/(1-rho))*(`ssalltc' + gamma*`myinicd') - lam`t'*`at'[1,nssrhs+T+2+`t']) `if'; qui replace mm`t'=`m`t''; * noisily scalar di rho; * noisily scalar di gamma; local t=`t'+1; }; end; ***************************************************************************; capture program drop gmm_dselect02; // Main equation, selection ignored; program gmm_dselect02; version 11; syntax varlist if, at(name) mysslhs(varlist) myinicd(varlist) myssrhs(varlist); local t=1; while `t'<=T {; local sslhs`t': word `t' of `mysslhs'; // second-step LHS; local t=`t'+1; }; // Obtain parameters on 2nd step time-constant (x1,...xt,z1,...,zt-educ) variables // and compute the total contribution of these variables to be // multiplied by rho function later, when obtaining moment conditions; tempvar ssalltc; qui gen double `ssalltc'=0 `if'; local k=0; foreach var of varlist `myssrhs' {; local k=`k'+1; qui replace `ssalltc'= `ssalltc' + `var'*`at'[1,`k'] `if'; }; scalar nsstcrhs=`k'; // Excludes the intercept; * noisily scalar di nsstcrhs; * noisily di "2nd step cum eff. of tc vars, ssalltc"; * noisily sum `ssalltc'; // Obtain parameters on second-step time-varying (x) variables // (following time-constant second-step parameters); foreach var of varlist $sstvvars {; local k=`k'+1; scalar beta`var'=`at'[1,`k']; * noisily scalar di beta`var'; }; scalar nssrhs=`k'; // Excludes a constant and initial condition; * noisily scalar di nssrhs; // Obtain parameter eta11...eta1T (following 2nd step x & z variables); local t=1; while `t'<=T {; scalar eta1`t'=`at'[1,nssrhs+`t']; local t=`t'+1; }; // Obtain parameter gamma (following 2nd step x & z variables and eta11...eta1T); scalar gamma=`at'[1,nssrhs+T+1]; // Obtain parameter rho (following 2nd step x & z variables, eta11...eta1T and gamma); scalar rho=`at'[1,nssrhs+T+2]; // Generate starting values for time-varying variables // to be used later for computing moment conditions at the 2nd step; foreach var of varlist $sstvvars {; // Accum. second-step time-varying X; tempvar cum`var'0; qui gen double `cum`var'0'=0 `if'; * noisily di "cumulative effect of tv variable" `var' ", t=0"; * noisily sum `cum`var'0'; }; tempvar cumcons0; qui gen double `cumcons0'=0 `if'; //Start the cycle over time - compute 2nd step moments for each t; local t=1; while `t'<=T {; * noisily di "Time period t = " `t'; // Add up all cumulative contributions of second-step time-varying X // at time t; local t_lag=`t'-1; tempvar ssalltv`t'; qui gen double `ssalltv`t''=0 `if'; foreach var of varlist $sstvvars {; tempvar cum`var'`t'; qui gen double `cum`var'`t''=rho*`cum`var'`t_lag''+`var'`t' `if'; qui replace `ssalltv`t''=`ssalltv`t''+`cum`var'`t''*beta`var' `if'; }; // Add the cumulative effect of time-specific intercept tempvar cumcons`t'; qui gen double `cumcons`t''=rho*`cumcons`t_lag''+eta1`t' `if'; qui replace `ssalltv`t''=`ssalltv`t''+`cumcons`t'' `if'; * noisily di "cum. contribution of all 2nd step tv X at time " `t'; * noisily sum `ssalltv`t''; // Generate 2-step moment condition for time period t; local m`t': word `t' of `varlist'; qui replace `m`t''= s`t'*(`sslhs`t'' - (rho^`t')*`myinicd' - `ssalltv`t'' - ((1-rho^`t')/(1-rho))*(`ssalltc' + gamma*`myinicd')) `if'; * noisily scalar di rho; * noisily scalar di gamma; local t=`t'+1; }; end; **********************************************************************; // Do GMM using selected sample only; gen cons=1; gen age=.; gen agesq=.; global sstvvars age agesq; ************************************************************************* * DO ONE-STEP GMM AND USE RESULTS TO OBTAIN THE OPTIMAL WEIGHTING MATRIX *************************************************************************; gmm gmm_dselect, nequations(12) nparameters(41) mysslhs(lnw*) myinicd(ini_lnw) myssrhs(educ children*) instruments(1: ini_lnw educ age1 agesq1 children* lam1) instruments(2: ini_lnw educ age2 agesq2 children* lam2) instruments(3: ini_lnw educ age3 agesq3 children* lam3) instruments(4: ini_lnw educ age4 agesq4 children* lam4) instruments(5: ini_lnw educ age5 agesq5 children* lam5) instruments(6: ini_lnw educ age6 agesq6 children* lam6) instruments(7: ini_lnw educ age7 agesq7 children* lam7) instruments(8: ini_lnw educ age8 agesq8 children* lam8) instruments(9: ini_lnw educ age9 agesq9 children* lam9) instruments(10: ini_lnw educ age10 agesq10 children* lam10) instruments(11: ini_lnw educ age11 agesq11 children* lam11) instruments(12: ini_lnw educ age12 agesq12 children* lam12) winitial(unadjusted, independent) onestep nocommonesample; *** Use onestep if no efficiency gain OR to reduce time requirements; mat b=e(b); mat S=e(S); mat V=e(V)'; mat b1=b[1,30..41]; mat V1=V[30..41,30..41]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2=chi2[1,1]; scalar pv_ch2=chi2tail(12,ch2); scalar di ch2; scalar di pv_ch2; scalar xi_educ=_b[/b1]; scalar b_age=_b[/b14]; scalar b_agesq=_b[/b15]; scalar gamma=_b[/b28]; scalar rho=_b[/b29]; local i=2; foreach var of varlist children* {; scalar xi_`var'=_b[/b`i']; scalar di xi_`var'; local i=`i'+1; }; local i=16; local t=1; while `i'<=27 {; scalar eta`t'=_b[/b`i']; local i=`i'+1; local t=`t'+1; }; local i=30; local t=1; while `i'<=41 {; scalar vp`t'=_b[/b`i']; local i=`i'+1; local t=`t'+1; }; *****************************************************************; * OBTAINING DERIVATIVES *****************************************************************; * DERIVATIVE WRT COEf. ON AGE, AGESQ; foreach var of varlist age agesq {; gen dv_b_`var'1=-s1*`var'1; }; local t=2; while `t'<=T {; local t_lag=`t'-1; foreach var of varlist age agesq {; gen dv_b_`var'`t'=-s`t'*(`var'`t'+rho*dv_b_`var'`t_lag'); }; local t=`t'+1; }; * DERIVATIVE WRT COEf. ON EDUC; local t=1; while `t'<=T {; gen dv_xi_educ`t'=-s`t'*(1-rho^`t')/(1-rho)*educ; local t=`t'+1; }; * DERIVATIVE WRT COEF. ON CHILDREN1-12; local t=1; while `t'<=T {; foreach var of varlist children* {; gen dv_xi_`var'_`t'=-s`t'*(1-rho^`t')/(1-rho)*`var'; }; local t=`t'+1; }; * DERIVATIVE WRT TIME-SPECIFIC INTERCEPT; local t=1; while `t'<=T {; local j=1; while `j'<=T {; gen dv_eta`t'_`j'=0; local j=`j'+1; }; local t=`t'+1; }; local t=1; while `t'<=T {; local j=1; while `j'<=`t' {; replace dv_eta`j'_`t'=-s`t'*rho^(`t'-`j'); local j=`j'+1; }; local t=`t'+1; }; * DERIVATIVE WRT GAMMA; local t=1; while `t'<=T {; gen dv_gamma`t'=-s`t'*(1-rho^`t')/(1-rho)*ini_lnw; local t=`t'+1; }; * DERIVATIVE WRT RHO; gen tv_dv_rho_1=0; local t=2; while `t'<=T {; gen tv_dv_rho_`t'=0; local j=1; local t_lag=`t'-1; while `j'<=`t_lag' {; local tj=`t'-`j'; replace tv_dv_rho_`t'=-s`t'*(tv_dv_rho_`t'+ `j'*rho^(`j'-1)*(age`tj'*b_age + agesq`tj'*b_agesq + eta`tj')); local j=`j'+1; }; local t=`t'+1; }; gen double all_children=0; foreach var of varlist children* {; replace all_children = all_children + xi_`var'*`var'; }; sum all_children; local t=1; while `t'<=T {; gen double dv_rho_`t'=-s`t'*[ `t'*rho^(`t'-1)*ini_lnw + tv_dv_rho_`t' + ( (1-rho^`t')/(1-rho)^2 - `t'*rho^(`t'-1)/(1-rho) )*(educ*xi_educ + all_children + gamma*ini_lnw) ]; local t=`t'+1; }; * DERIVATIVE WRT PSI (COEf. ON LAMBDA); local t=1; while `t'<=T {; local j=1; while `j'<=T {; gen dv_psi`t'_`j'=0; local j=`j'+1; }; local t=`t'+1; }; local t=1; while `t'<=T {; replace dv_psi`t'_`t'=-s`t'*lam`t'; local t=`t'+1; }; *****************************************************************; * COMPUTING MATRIX G = W' * *****************************************************************; * 18 = # IVs (at 2nd step) in period t, incl. cons - same in each * 41 = # parameters to estimate at 2nd step *****************************************************************; mat G=J(18*T,41,0); local t=1; while `t'<=T {; local W`t' ini_lnw educ age`t' agesq`t' children* lam`t' cons; local dv_g`t' dv_xi_educ`t' dv_xi_children1_`t' dv_xi_children2_`t' dv_xi_children3_`t' dv_xi_children4_`t' dv_xi_children5_`t' dv_xi_children6_`t' dv_xi_children7_`t' dv_xi_children8_`t' dv_xi_children9_`t' dv_xi_children10_`t' dv_xi_children11_`t' dv_xi_children12_`t' dv_b_age`t' dv_b_agesq`t' dv_eta1_`t' dv_eta2_`t' dv_eta3_`t' dv_eta4_`t' dv_eta5_`t' dv_eta6_`t' dv_eta7_`t' dv_eta8_`t' dv_eta9_`t' dv_eta10_`t' dv_eta11_`t' dv_eta12_`t' dv_gamma`t' dv_rho_`t' dv_psi1_`t' dv_psi2_`t' dv_psi3_`t' dv_psi4_`t' dv_psi5_`t' dv_psi6_`t' dv_psi7_`t' dv_psi8_`t' dv_psi9_`t' dv_psi10_`t' dv_psi11_`t' dv_psi12_`t'; mat accum TEMP= `W`t'' `dv_g`t'', nocons; mat G[18*(`t'-1)+1,1]=TEMP[1..18,19...]; local t=`t'+1; }; *****************************************************************; * OBTAINING CORRECTED STD. ERRORS FOR ONE-STEP GMM *****************************************************************; *****************************************************************; * varli2, varli3 AND varli4 ARE EMPTY FOR NOW AND WILL BE USED LATER *****************************************************************; local t=1; while `t'<=T {; local varli2`t'; local varli3`t'; local varli4`t'; local t=`t'+1; }; mat H=I(17*12); /* 17 = # vars in probit, 12=T */ local t=1; while `t'<=T {; qui probit s`t' age`t' agesq`t' educ ini_lnw children*; qui predict xb, xb; mat H`t'=e(V); mat H[17*(`t'-1)+1,17*(`t'-1)+1]=H`t'; qui gen tempvar1=normalden(xb)/normal(xb)*s`t'; assert tempvar1==lam`t'; qui replace tempvar1=-normalden(xb)/(1-normal(xb)) if s`t'==0; foreach var of varlist age`t' agesq`t' educ ini_lnw children* cons {; qui gen double g`var'_`t' = -s`t'*`var'*vp`t'*lam`t'*(lam`t'+xb); qui gen double q`var'_`t' = tempvar1*`var'; local varli3`t' `varli3`t'' g`var'_`t'; local varli4`t' `varli4`t'' q`var'_`t'; }; drop xb tempvar*; local t=`t'+1; }; *****************************************************************; * DEFINE varli2 - THE LIST OF INTERACTION TERMS * (Z*) *****************************************************************; local t=1; while `t'<=T {; foreach var of varlist ini_lnw educ age`t' agesq`t' children* lam`t' cons {; qui gen double eh`var'_`t'=`var'*mm`t'; // NOTE: mm`t'=0 if s`t'=0; local varli2`t' `varli2`t'' eh`var'_`t'; }; local t=`t'+1; }; *****************************************************************; * FOR F BELOW, 18 = # IV vars at 2nd step, 17 = # vars at 1st step *****************************************************************; matrix F=J(T*18,T*17,0); local t=1; while `t'<=T {; mat accum TEMP=ini_lnw educ age`t' agesq`t' children* lam`t' cons `varli3`t'', nocons; matrix F`t'=TEMP[1..18,19...]; matrix F[18*(`t'-1)+1,17*(`t'-1)+1]=F`t'; local t=`t'+1; }; matrix F=F/`nobs'; local varli2T; local varli4T; local t=1; while `t'<=T {; local varli2T `varli2T' `varli2`t''; local varli4T `varli4T' `varli4`t''; local t=`t'+1; }; mat accum WggW=`varli2T', nocons; mat accum dd=`varli4T', nocons; mat accum TEMP=`varli2T' `varli4T', nocons; mat Wgd=TEMP[1..18*T,18*T+1...]; mat CS=(WggW-F*H*Wgd'-(F*H*Wgd')'+F*H*dd*H*F')/`nobs'; mat WggW_N=WggW/`nobs'; *mat list CS; *mat list S; *mat list WggW_N; *****************************************************************; * CORRECTED STD. ERRORS FOR ONE-STEP GMM WITH ESTIMATED LAMBDA *****************************************************************; mat V_NLS=syminv(G'*G)*G'*CS*G*syminv(G'*G)*`nobs'; scalar se_educ=sqrt(V_NLS[1,1]); scalar se_age=sqrt(V_NLS[14,14]); scalar se_agesq=sqrt(V_NLS[15,15]); scalar se_rho=sqrt(V_NLS[29,29]); scalar di se_educ; scalar di se_age; scalar di se_agesq; scalar di se_rho; mat V=V_NLS; mat b1=b[1,30..41]; mat V1=V[30..41,30..41]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2=chi2[1,1]; scalar pv_ch2=chi2tail(12,ch2); scalar di ch2; scalar di pv_ch2; ************************************************************************* * DO TWO-STEP GMM USING THE OPTIMAL WEIGHTING MATRIX *************************************************************************; mat OW=syminv(CS); gmm gmm_dselect, nequations(12) nparameters(41) mysslhs(lnw*) myinicd(ini_lnw) myssrhs(educ children*) instruments(1: ini_lnw educ age1 agesq1 children* lam1) instruments(2: ini_lnw educ age2 agesq2 children* lam2) instruments(3: ini_lnw educ age3 agesq3 children* lam3) instruments(4: ini_lnw educ age4 agesq4 children* lam4) instruments(5: ini_lnw educ age5 agesq5 children* lam5) instruments(6: ini_lnw educ age6 agesq6 children* lam6) instruments(7: ini_lnw educ age7 agesq7 children* lam7) instruments(8: ini_lnw educ age8 agesq8 children* lam8) instruments(9: ini_lnw educ age9 agesq9 children* lam9) instruments(10: ini_lnw educ age10 agesq10 children* lam10) instruments(11: ini_lnw educ age11 agesq11 children* lam11) instruments(12: ini_lnw educ age12 agesq12 children* lam12) winitial(OW) onestep vce(unadjusted) from(b) nocommonesample; scalar b_educ_GMM_lam=_b[/b1]; scalar b_age_GMM_lam=_b[/b14]; scalar b_agesq_GMM_lam=_b[/b15]; scalar rho_GMM_lam=_b[/b29]; scalar se_educ_GMM_lam=_se[/b1]; scalar se_age_GMM_lam=_se[/b14]; scalar se_agesq_GMM_lam=_se[/b15]; scalar se_rho_GMM_lam=_se[/b29]; mat b=e(b); mat V=e(V)'; mat b1=b[1,30..41]; mat V1=V[30..41,30..41]; mat chi2=b1*inv(V1)*b1'; mat list chi2; scalar ch2_GMM_lam=chi2[1,1]; scalar pv_ch2_GMM_lam=chi2tail(12,ch2_GMM_lam); scalar di ch2_GMM_lam; scalar di pv_ch2_GMM_lam; ***************************************************************** * SUMMARIZE ALL RESULTS *****************************************************************; * FD-IV W/O LAMBDA; scalar di b_agesq_FDIV_nlam; scalar di rho_FDIV_nlam; scalar di se_agesq_FDIV_nlam; scalar di se_rho_FDIV_nlam; *** FD-IV WITH LAMBDA; scalar di b_agesq_FDIV_lam; scalar di rho_FDIV_lam; scalar di se_agesq_FDIV_lam; scalar di se_rho_FDIV_lam; scalar di ch2_FDIV; scalar di pv_ch2_FDIV; *** NLS W/O LAMBDA; scalar di b_educ_NLS_nlam; scalar di b_age_NLS_nlam; scalar di b_agesq_NLS_nlam; scalar di rho_NLS_nlam; scalar di se_educ_NLS_nlam; scalar di se_age_NLS_nlam; scalar di se_agesq_NLS_nlam; scalar di se_rho_NLS_nlam; *** NLS WITH LAMBDA; scalar di b_educ_NLS_lam; scalar di b_age_NLS_lam; scalar di b_agesq_NLS_lam; scalar di rho_NLS_lam; scalar di se_educ_NLS_lam; scalar di se_age_NLS_lam; scalar di se_agesq_NLS_lam; scalar di se_rho_NLS_lam; * USE VAR NOT CORRECTED FOR 1st STEP ESTIMATION; scalar di ch2_NLS; scalar di pv_ch2_NLS; * USE VAR CORRECTED FOR 1st STEP ESTIMATION; scalar di ch2_NLS_lam; scalar di pv_ch2_NLS_lam; *** TWO-STEP GMM WITH LAMBDA; scalar di b_educ_GMM_lam; scalar di b_age_GMM_lam; scalar di b_agesq_GMM_lam; scalar di rho_GMM_lam; scalar di se_educ_GMM_lam; scalar di se_age_GMM_lam; scalar di se_agesq_GMM_lam; scalar di se_rho_GMM_lam; scalar di ch2_GMM_lam; scalar di pv_ch2_GMM_lam; log close;