olsnw.m (2563B)
1 % OLSNW.m -- Least Squares procedure -- with Newey-West standard errors 2 % 2/27/2001 -- merged ols.m and lsnw.m 3 % [u,N,K,stdest,beta,vcv,robvcv] = ols(y,x,title,depv,indv,prevest); 4 % prevest = degrees of freedom correction = # of parameters previously est. 5 % 11/93 - Deletes missing values if encountered. 6 7 function [u,N,K,stdest,beta,vcv,robvcv] = olsnw(y,x,title,depv,indv,prevest); 8 9 if x(:,1) ~= 1; disp 'No constant term in regression'; end; 10 if exist('prevest')~=1; prevest=0; end; 11 if title == 0; title = 'Ordinary Least Squares: Newey-West S.E.s'; end; 12 if depv == 0; depv = ' '; end; 13 14 [N K] = size(x); 15 if indv == 0; indv=vdummy('x',K); end; 16 17 % Check for missing values 18 data=packr([y x]); 19 if size(data,1)~=N; disp 'Missing values encountered'; 20 y=data(:,1); 21 x=data(:,2:K+1); 22 [N K] = size(x); 23 end; 24 25 xxinv = inv(x'*x); 26 beta = xxinv*x'*y; 27 u = y-x*beta; 28 dof = N-K-prevest; 29 sigma2=u'*u/dof; 30 stdest=sqrt(sigma2); 31 vcv =sigma2*xxinv; 32 se =sqrt(diag(vcv)); 33 tstat = beta./se; 34 ybar = 1/N*(ones(N,1)'*y); 35 rsq = 1-(u'*u)/(y'*y - N*ybar^2); 36 rbar = 1-sigma2/((y'*y - N*ybar^2)/(N-1)); 37 38 % Compute Newey-West standard errors using lsnw 39 [bb,trob,told,ssigma2,robvcv] = lsnw(y,x); 40 roberr = sqrt(diag(robvcv)); 41 42 disp '============================================================================='; 43 disp ' '; 44 disp ' '; 45 disp ' ------------------------------------------------------------'; 46 fprintf([' ', title, '\n']); 47 disp ' ------------------------------------------------------------'; 48 disp ' '; 49 fprintf([' NOBS: ', num2str(N),' Dependent Variable: ',depv,'\n']); 50 fprintf([' RHS Vars: ', num2str(K),'\n']); 51 fprintf([' D of F: ', num2str(dof), '\n']); 52 disp ' '; 53 fprintf([' R-Squared: ', num2str(rsq)]); 54 fprintf([' S.E.E.: ', num2str(stdest), '\n']); 55 fprintf([' RBar^2: ', num2str(rbar)]); 56 fprintf([' Residual SS: ', num2str(u'*u), '\n']); 57 disp ' '; 58 disp ' Standard NewWes NewWes'; 59 disp 'Variable Beta Error t-stat Error t-stat'; 60 disp '-------- ---- -------- ------ ------ ------'; 61 disp ' '; 62 fmt='%16.6f %12.6f %8.2f %12.6f %8.2f'; 63 results=[beta se tstat roberr trob]; 64 for i=1:K; 65 fprintf(indv(i,:)); 66 fprintf(1,fmt,results(i,:)); 67 fprintf('\n'); 68 end 69 disp ' '; 70 disp ' '; 71 disp '============================================================================='; 72 %end; % if PRINT==0 73 74 K=K+prevest; % Return Correct # of Estimated Pars