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