/* example of testing in SUR system consider a two equation system for simplicity, with two regressors in each equation */ output file = sur_example.out reset; /* sample size is n */ n_vals = {20,50,400}; jjj = 1; do until jjj gt 3; n = n_vals[jjj]; x1_1 = rndu(n,1); x2_1 = 2 + .3*x1_1 + rndu(n,1); x1_2 = rndu(n,1); x2_2 = .6*x1_2 + rndu(n,1); x1 = ones(n,1) ~ x1_1 ~ x2_1; x2 = ones(n,1) ~ x1_2 ~ x2_2; x = zeros(2*n,6); x[1:n,1:3] = x1; s2 = n+1; f2 = 2*n; x[s2:f2,4:6] = x2; /* first assume that disturbances are normal with the usual SUR structure - generate R replications */ r = 2000; /* this section constructs normally distributed disturbances */ print "**********************************"; print "THIS RUN USES NORMAL DISTURBANCES"; PRINT "**********************************"; @true sigma matrix@ sigma = {4 3, 3 6}; c = chol(sigma); eps_mat = zeros(2*n,r); i = 1; do until i gt r; jj = 1; do until jj gt n; eps_mat[jj (jj+n),i] = c'rndn(2,1); jj = jj+1; endo; i = i+1; endo; eps_1 = eps_mat[1:n,.]; eps_2 = eps_mat[s2:f2,.]; /* end of construction of normal disturbances */ /* this section constructs non normal disturbances we try to match the first two moments of the normal disturbances as closely as possible */ /*********** print "**********************************"; print "THIS RUN USES NON-NORMAL DISTURBANCES"; print " BASED ON UNIFORM "; print "**********************************"; eps_1 = sqrt(48)*(rndu(n,r)-.5); eps_2 = .75*eps_1 + 6.7*(rndu(n,r)-.5); /* "true" sigma in this case */ sigma = {4 2.979, 2.979 5.961}; ************/ /* this ends the generation of uniform-based disturbances */ /* actual parameter values */ beta_1 = {.4,3,0}; beta_2 = {.7,0,1}; /* generate data */ y1_mat = x1*beta_1 + eps_1; y2_mat = x2*beta_2 + eps_2; /* compute ols system estimator */ beta_ols = zeros(6,r); i = 1; do until i gt r; y = y1_mat[.,i] | y2_mat[.,i]; beta_ols[.,i] = invpd(x'x)*x'y; i = i+1; endo; /* compute mean of ols estimates and their standard errors */ /* estimate sigma matrix in SUR setup */ sigma_hat = zeros(3,r); i = 1; do until i gt r; y = y1_mat[.,i] | y2_mat[.,i]; res = y - x*beta_ols[.,i]; sigma_hat[1,i] = res[1:n]'res[1:n]/(n-3); sigma_hat[2,i] = res[1:n]'res[s2:f2]/(n-3); sigma_hat[3,i] = res[s2:f2]'res[s2:f2]/(n-3); i = i+1; endo; /* tests of restriction on parameters */ cap_r = {0 0 1 0 0 0, 0 0 0 0 1 0}; small_r = {0,0}; /* assume that sigma is known */ test_known_sigma = zeros(r,1); i = 1; /* matrix of quadratic form in test statistic */ qform_known = invpd(cap_r*(invpd(x'x)*x'(sigma .*. eye(n))*x*invpd(x'x))*cap_r'); do until i gt r; test_known_sigma[i] = (cap_r*beta_ols[.,i]-small_r)'qform_known*(cap_r*beta_ols[.,i]-small_r); i = i+1; endo; library pgraph; graphset; @histp(test_known_sigma,50);@ /* now compute test using OLS estimator with estimated sigma */ test_unknown_sigma = zeros(r,1); i = 1; do until i gt r; sigma_r = zeros(2,2); sigma_r[1,1] = sigma_hat[1,i]; sigma_r[1,2] = sigma_hat[2,i]; sigma_r[2,1] = sigma_r[1,2]; sigma_r[2,2] = sigma_hat[3,i]; qform_unknown_r = invpd(cap_r*(invpd(x'x)*x'(sigma_r .*. eye(n))*x*invpd(x'x))*cap_r'); test_unknown_sigma[i] = (cap_r*beta_ols[.,i]-small_r)'qform_unknown_r*(cap_r*beta_ols[.,i]-small_r); i = i+1; endo; /*********** GLS Part ***********/ beta_gls = zeros(6,r); i = 1; psi_inverse = invpd(sigma) .*. eye(n); do until i gt r; y = y1_mat[.,i] | y2_mat[.,i]; beta_gls[.,i] = invpd(x'psi_inverse*x)*x'psi_inverse*y; i = i+1; endo; /* assume that sigma is known */ test_known_sigma_gls = zeros(r,1); i = 1; /* matrix of quadratic form in test statistic */ qform_known_gls = invpd(cap_r*invpd(x'(invpd(sigma) .*. eye(n))*x)*cap_r'); do until i gt r; test_known_sigma_gls[i] = (cap_r*beta_gls[.,i]-small_r)'qform_known_gls*(cap_r*beta_gls[.,i]-small_r); i = i+1; endo; library pgraph; graphset; @histp(test_known_sigma_gls,50);@ /*** FEASIBLE GLS PART *****/ beta_fgls = zeros(6,r); /* in order to estimate sigma, use ols residuals from above, and the OLS-based estimates of sigma */ i = 1; do until i gt r; sigma_r = zeros(2,2); sigma_r[1,1] = sigma_hat[1,i]; sigma_r[1,2] = sigma_hat[2,i]; sigma_r[2,1] = sigma_r[1,2]; sigma_r[2,2] = sigma_hat[3,i]; y = y1_mat[.,i] | y2_mat[.,i]; beta_fgls[.,i] = invpd(x'(invpd(sigma_r) .*. eye(n))*x)*x'(invpd(sigma_r) .*. eye(n))*y; i = i+1; endo; /* now compute test using FGLS estimator with estimated sigma */ test_unknown_sigma_fgls = zeros(r,1); i = 1; do until i gt r; sigma_r = zeros(2,2); sigma_r[1,1] = sigma_hat[1,i]; sigma_r[1,2] = sigma_hat[2,i]; sigma_r[2,1] = sigma_r[1,2]; sigma_r[2,2] = sigma_hat[3,i]; qform_unknown_gls_r = invpd(cap_r*invpd(x'(invpd(sigma_r) .*. eye(n))*x)*cap_r'); test_unknown_sigma_fgls[i] = (cap_r*beta_fgls[.,i]-small_r)'qform_unknown_gls_r*(cap_r*beta_fgls[.,i]-small_r); i = i+1; endo; print ""; print "***************** Parameter Estimates ******************/"; print "SAMPLE SIZE = " n; print " estimates of covariance matrix "; print " true value mean st dev "; print ((sigma[1,1] | sigma[1,2] | sigma[2,2]) ~ meanc(sigma_hat') ~ stdc(sigma_hat')); print "************************"; print " TRUE VALUES OLS mean OLS st dev GLS mean GLS st dev FGLS mean FGLS st dev"; print ((beta_1 | beta_2) ~ meanc(beta_ols') ~ stdc(beta_ols') ~ meanc(beta_gls') ~ stdc(beta_gls') ~ meanc(beta_fgls') ~ stdc(beta_fgls')); /* test statistic distributions compute 1 - cdf(evaluation point) */ epoint = seqa(1,1,10); prob_mat = zeros(10,5); kk = 1; do until kk gt 10; prob_mat[kk,1] = cdfchic(epoint[kk],2); prob_mat[kk,2] = meanc(test_known_sigma .> epoint[kk]); prob_mat[kk,3] = meanc(test_unknown_sigma .> epoint[kk]); prob_mat[kk,4] = meanc(test_known_sigma_gls .> epoint[kk]); prob_mat[kk,5] = meanc(test_unknown_sigma_fgls .> epoint[kk]); kk = kk+1; endo; /* print test statistic distributions */ print ""; print "************** Test statistics *****************"; print " eval point true prob OLS sigma known OLS sigma est GLS sigma known FGLS sigma est "; print (epoint ~ prob_mat); jjj = jjj+1; endo; end;