ivmiss.m (3390B)
1 % ivmiss [u,N,K,stdest,b,vcv,sOID]=ivmiss(y,x,w,title,depv,indv,instv) 2 % Instrumental Variables Estimation (Works with Panel data also) 3 % 4 % y = dependent variable Nx1 5 % x = rhs variables NxK 6 % w = instruments NxJ J>=K 7 % 8 % ivmiss drops missing data from [y x w] 7/9/12 9 10 function [u,N,K,stdest,beta,vcv,sOID]=ivmiss(y,x,w,title,depv,indv,instv); 11 12 sOID=[]; 13 14 % Drop missing data first 15 keep=~any(isnan([y x w]'))'; 16 y=y(keep); 17 x=x(keep,:); 18 w=w(keep,:); 19 20 Pw = w*inv(w'*w)*w'; % The projection matrix for w 21 xxinv = inv(x'*Pw*x); 22 beta = xxinv*x'*Pw*y; 23 u = y-x*beta; % Use the real x's here! 24 25 % Now, just copy the remainder from ols.m 26 27 if x(:,1) ~= 1; disp 'No constant term in regression'; end; 28 if exist('prevest')~=1; prevest=0; end; 29 if title == 0; title = 'Instrumental Variables Estimation'; end; 30 if depv == 0; depv = ' '; end; 31 32 [N K] = size(x); 33 if indv == 0; indv=vdummy('x',K); end; 34 35 36 dof = N-K-prevest; 37 sigma2=u'*u/dof; 38 stdest=sqrt(sigma2); 39 vcv =sigma2*xxinv; 40 se =sqrt(diag(vcv)); 41 tstat = beta./se; 42 ybar = 1/N*(ones(N,1)'*y); 43 rsq = 1-(u'*u)/(y'*y - N*ybar^2); 44 rbar = 1-sigma2/((y'*y - N*ybar^2)/(N-1)); 45 46 % Compute White robust standard errors 47 48 robvcv=zeros(K,K); 49 xhat=Pw*x; 50 for i=1:N 51 robvcv = robvcv + u(i)^2*xhat(i,:)'*xhat(i,:); 52 end % i 53 robvcv = (N/(dof))*xxinv*robvcv*xxinv; 54 roberr = sqrt(diag(robvcv)); 55 trob = beta./roberr; 56 57 58 % Test of OID Restrictions (See Econometrics Notes, Newey Section (end) back 59 % of page 5: Based on quadratic form for W'e/sqrt(T). 60 61 dofOID=size(w,2)-size(x,2); 62 if dofOID>0; 63 sOID=(w'*u)'*inv(sigma2*w'*w)*(w'*u); 64 pval=1-chi2cdf(sOID,dofOID); 65 end; 66 67 % Hausman-Wu test versus OLS; Problem: often VV is not invertible (pd) 68 %bols = inv(x'*x)*x'*y; 69 %q=bols-beta; 70 %VV=sigma2*(xxinv - inv(x'*x)); 71 %det(VV) 72 %sHW=q'*inv(VV)*q; 73 %pHW=1-chi2cdf(sHW,K); 74 75 disp '============================================================================='; 76 disp ' '; 77 disp ' '; 78 disp ' ------------------------------------------------------------'; 79 fprintf([' ', title, '\n']); 80 disp ' ------------------------------------------------------------'; 81 disp ' '; 82 fprintf([' NOBS: ', num2str(N),' Dependent Variable: ',depv,'\n']); 83 fprintf([' RHS Vars: ', num2str(K),'\n']); 84 fprintf([' D of F: ', num2str(dof), '\n']); 85 disp ' '; 86 fprintf([' R-Squared: ', num2str(rsq)]); 87 fprintf([' S.E.E.: ', num2str(stdest), '\n']); 88 fprintf([' RBar^2: ', num2str(rbar)]); 89 fprintf([' Residual SS: ', num2str(u'*u), '\n']); 90 disp ' '; 91 if dofOID>0; 92 fprintf(' OverID Test -- Statistic: %6.2f DoF: %4.0f Pval: %6.3f\n',[sOID dofOID pval]); 93 end; 94 %fprintf(' Hausman-Wu Test -- Stat: %6.2f DoF: %4.0f Pval: %6.3f\n',[sHW K pHW]); 95 disp ' '; 96 disp ' Standard Robust Robust'; 97 disp 'Variable Beta Error t-stat Error t-stat'; 98 disp '-------- ---- -------- ------ ------ ------'; 99 disp ' '; 100 fmt='%16.6f %12.6f %8.2f %12.6f %8.2f'; 101 results=[beta se tstat roberr trob]; 102 for i=1:K; 103 fprintf(indv(i,:)); 104 fprintf(1,fmt,results(i,:)); 105 fprintf('\n'); 106 end 107 disp ' '; 108 fprintf('Instruments: '); 109 say(instv); % List the instruments 110 disp ' '; 111 disp '============================================================================='; 112 %end; % if PRINT==0