Wald test in Stata

* Note that x = 1993 and y = 1998.

cap drop Insample
mark Insample if logpay93 !=. & logpay98 != .
summ Insample
scalar N = r(sum)

* 1st entry of covariance matrix
corr logpay93 logpay98, cov
scalar var_mu93_minus_mu98 = (r(Var_1)*(r(N) - 1)/r(N) + r(Var_2)*(r(N)-1)/r(N) - 2*r(cov_12)*(r(N)-1)/r(N))/N
// Stata computes the sample estimates, but we are working with MLEs

* 2nd and 3rd entries (cov_mus_vars)
foreach y in 93 98 {
summ logpay`y', detail
scalar mu_`y' = r(mean)
scalar var_`y' = r(Var)*(r(N) - 1)/r(N)
gen logpay`y'deviation2 = (logpay`y' - mu_`y')^2

correlate logpay`y' logpay`y'deviation2, covariance
scalar cov_`y'_`y'deviation2 = r(cov_12)*(r(N) - 1)/r(N)
scalar var_`y'deviation2 = r(Var_2)*(r(N)-1)/r(N)  // comes in handy for the final entry
display cov_`y'_`y'deviation2
}

corr logpay93 logpay98deviation2, cov
scalar cov_93_98deviation2 = r(cov_12)*(r(N) - 1)/r(N)
display cov_93_98deviation2

corr logpay98 logpay93deviation2, cov
scalar cov_98_93deviation2 = r(cov_12)*(r(N) - 1)/r(N)
display cov_98_93deviation2

scalar cov_mus_vars = ( cov_93_93deviation2 - cov_93_98deviation2 - cov_98_93deviation2 + cov_98_98deviation2 )/N
display  cov_mus_vars


* 4th entry
correlate logpay93deviation2  logpay98deviation2, cov
scalar cov_dev2s = r(cov_12)*(r(N) - 1)/r(N)

scalar var_vars = ( var_93deviation2 + var_98deviation2 - cov_dev2s )/N

di var_mu93_minus_mu98
cap drop Cov_theta
matrix define Cov_theta = ( var_mu93_minus_mu98 , cov_mus_vars \  cov_mus_vars, var_vars )
matrix list Cov_theta, nohalf // note that only the lower triangle will be printed if "nohalf" is not specified, as this matrix is symmetric

matrix inv_Cov_theta = invsym(Cov_theta)
matrix list inv_Cov_theta

matrix define theta = ( mu_93 - mu_98 \ var_93 - var_98 )
mat thetarow = theta'

matrix Wald_joint_matrix = thetarow*inv_Cov_theta*theta
scalar Wald_joint = Wald_joint_matrix[1, 1]
di Wald_joint // Wald statistic
di chi2tail(2, Wald_joint) // p-value