panelreg.m (8019B)
1 % panelreg.m function [u,NN,K,stdest,beta,vcv]= 2 % panelreg(y,x,Method,yname,xnames,nnames); 3 % This is the FUNCTION for the panel regressions. 4 % Run the panel regressions in levels including fixed effects: 5 % 6 % y = a_i +b_i*t + Xbeta + epsilon 7 % 8 % and then we can do some hypothesis testing. 9 % 10 % Method is a 2x1 vector Method=[FE Time] 11 % FE=0 ==> No Fixed Effects 12 % FE=1 ==> Fixed Effects 13 % FE=2 ==> F.E. and All of the X coefficients vary. 14 % FE=3 ==> Between Estimator 15 % FE=4 ==> Random Effects--GLS (Theta differenced) 16 % Time=0 ==> No time trends 17 % Time=1 ==> Single time trend 18 % Time=2 ==> Country Specific Time trends (automatically FE) 19 % 20 % Note: If Method is scalar, Time=0 is assumed. Also, the between 21 % (FE=3) and the GLS (FE=4) assume no time trends... 22 23 24 function [u,NN,K,stdest,beta,vcv]=... 25 panelreg(y,x,Method,yname,xnames,nnames); 26 27 if length(Method)==1; 28 Time=0; 29 else; 30 Time=Method(2); 31 end; 32 FE=Method(1); 33 34 % So now we've got our data and the years, feed this into panlag. 35 % Take out means (mean0=1) so that we are running the fixed effects 36 % regression. Then we have to adjust the degrees of freedom. 37 % Notice: PanShape returns the Balanced Panel only. 38 39 [T0,N0]=size(y); % Years x Countries 40 [Tx,NX]=size(x); % Years x Countries*Variables 41 numv=NX/N0; % Number of Variables 42 X=[]; 43 44 % Random Effects: Greene, p. 490 45 % Run the Within and Between regressions and use the output to construct 46 % estimates of theta. Then call Panshape with Theta to 47 % theta-difference the data and then run the regression. 48 49 if FE==4; % Random Effects/GLS 50 51 % Pooled; 52 Meth=0; [uP,NP,KP,sigP,bP,vcvP]=panelreg(y,x,Meth,yname,xnames,nnames); 53 54 % Within; 55 Meth=1; [uW,NW,KW,sigW,bW,vcvW]=panelreg(y,x,Meth,yname,xnames,nnames); 56 Ftest(uP,uW,KW-KP,NW-KW,'Test H0: Fixed Effects are not different'); 57 58 % Between; 59 Meth=3; [uB,NB,KB,sigB,bB,vcvB]=panelreg(y,x,Meth,yname,xnames,nnames); 60 61 T=NW/NB; 62 sig2e=sigW^2; 63 sig2u=sigB^2-sig2e/T; 64 fprintf('Estimated Sig2u: %12.8f\n',sig2u); 65 fprintf('Estimated T*Sig2b: %12.8f\n',T*sigB^2); 66 fprintf('Estimated Sig2e: %12.8f\n',sig2e); 67 theta=1-sigW/sqrt(T*sig2u+sig2e); 68 fprintf('Estimated Theta: %7.4f\n',theta); 69 if theta<0; 70 disp 'NEGATIVE THETA using Sigma2(Between) ==> Trying Pooled'; 71 den=sum((sum(reshape(uP,T,NB)).^2)')/(NB-KP)/T; 72 fprintf('Estimated T*Sig2mu: %12.8f\n',den); 73 fprintf('Estimated Sig2e: %12.8f\n',sig2e); 74 theta=1-sigW/sqrt(den); 75 fprintf('Estimated ThetaPooled: %7.4f\n',theta); 76 end; 77 fetime=4; 78 else; 79 fetime=(FE>0)+(Time==2); % Param for panshape 1=demean 2=detrend 80 theta=[]; 81 end; 82 83 if FE==3; fetime=3; end; % Panshape =3 ==> between. 84 [y keeps]=panshape(y,fetime,theta); 85 for i=1:numv; 86 [xx Nx]=panshape(x(:,(i-1)*N0+(1:N0)),fetime,theta); 87 X=[X xx]; 88 keeps=keeps.*Nx; 89 end; 90 NN=sum(keeps); 91 prevest=0; 92 93 data=packr([y X]); 94 95 if data==[]; 96 disp 'Cannot run the regression -- data matrix is empty'; 97 else; 98 y=data(:,1); 99 data(:,1)=[]; 100 NT=length(y); 101 102 disp ' '; 103 fprintf('Observations Kept in the Sample: %5.0f\n',sum(keeps)); 104 say(nnames(keeps,:)); 105 disp ' '; 106 fprintf('Eliminated for NaN Reasons: %5.0f\n',sum(~keeps)); 107 if sum(~keeps)~=0; say(nnames(~keeps,:)); end; 108 disp ' '; 109 110 if FE==0 & Time==0; 111 rhs=[ones(NT,1) data]; 112 indv=['Constant'; padspace(xnames,8)]; 113 tle='Panel Estimation--Pooled Regression'; 114 elseif FE==3; 115 rhs=[ones(NT,1) data]; 116 indv=['Constant'; padspace(xnames,8)]; 117 tle='Panel Estimation--Between Estimate'; 118 elseif FE==1 & Time==0; 119 rhs=data; 120 indv=xnames; 121 tle='Panel Estimation--Fixed Effects'; 122 prevest=NN; 123 elseif FE==4; 124 rhs=[ones(NT,1)*(1-theta) data]; % Difference the Constant! 125 % rhs=[ones(NT,1) data]; % Difference the Constant! 126 indv=['Constant'; padspace(xnames,8)]; 127 tle='Panel Estimation--Random Effects GLS'; 128 prevest=NN-1; 129 elseif FE==0 & Time==1; 130 commont=kron(ones(NN,1),(1:(NT/NN))'); 131 rhs=[ones(NT,1) commont data]; 132 indv=['Constant'; 'Time '; padspace(xnames,8)]; 133 tle='Panel Estimation--No F.E.'; 134 elseif FE==1 & Time==1; 135 commont=kron(ones(NN,1),(1:(NT/NN))'); 136 rhs=[commont data]; 137 indv=['Time '; padspace(xnames,8)]; 138 tle='Panel Estimation -- F.E./CommonTrend'; 139 prevest=NN; 140 elseif FE==1 & Time==2; 141 rhs=data; 142 indv=xnames; 143 tle='Panel Estimation -- F.E./CountryTrends'; 144 prevest=2*NN; 145 else; 146 disp 'Not Yet Implemented!!!'; break; 147 end; 148 if exist('yname')==1; depv=yname; else; depv='Y '; end; 149 150 % Now, we would normally call OLS. However, instead, let's include OLS here 151 % explicitly. The motivation for this is that the standard errors for GLS 152 % have to be adjusted. 153 154 % [u,N,K,stdest,beta,vcv] =ols(y,rhs,tle,depv,indv,prevest); 155 x=rhs; 156 title=tle; 157 158 %%%%%%%%%%%%%%% Beg of OLS %%%%%%%%%%%%%%%% 159 if x(:,1) ~= 1; disp 'No constant term in regression'; end; 160 if exist('prevest')~=1; prevest=0; end; 161 if title == 0; title = 'Ordinary Least Squares'; end; 162 if depv == 0; depv = ' '; end; 163 164 [N K] = size(x); 165 if indv == 0; indv=vdummy('x',K); end; 166 167 % Check for missing values 168 data=packr([y x]); 169 170 171 if size(data,1)~=N; disp 'Missing values encountered'; 172 y=data(:,1); 173 x=data(:,2:K+1); 174 [N K] = size(x); 175 end; 176 177 xxinv = inv(x'*x); 178 beta = xxinv*x'*y; 179 u = y-x*beta; 180 dof = N-K-prevest; 181 sigma2=u'*u/dof; 182 stdest=sqrt(sigma2); 183 if FE~=4; 184 vcv =sigma2*xxinv; 185 else; 186 vcv =sigW^2*xxinv; 187 end; 188 se =sqrt(diag(vcv)); 189 tstat = beta./se; 190 ybar = 1/N*(ones(N,1)'*y); 191 rsq = 1-(u'*u)/(y'*y - N*ybar^2); 192 rbar = 1-sigma2/((y'*y - N*ybar^2)/(N-1)); 193 194 disp '============================================================================='; 195 disp ' '; 196 disp ' '; 197 disp ' ------------------------------------------------------------'; 198 fprintf([' ', title, '\n']); 199 disp ' ------------------------------------------------------------'; 200 disp ' '; 201 fprintf([' NOBS: ', num2str(N),' Dependent Variable: ',depv,'\n']); 202 fprintf([' RHS Vars: ', num2str(K),'\n']); 203 fprintf([' D of F: ', num2str(dof), '\n']); 204 disp ' '; 205 fprintf([' R-Squared: ', num2str(rsq)]); 206 fprintf([' S.E.E.: ', num2str(stdest), '\n']); 207 fprintf([' RBar^2: ', num2str(rbar)]); 208 fprintf([' Residual SS: ', num2str(u'*u), '\n']); 209 disp ' '; 210 disp ' Standard Robust Robust'; 211 disp 'Variable Beta Error t-stat Error t-stat'; 212 disp '-------- ---- -------- ------ ------ ------'; 213 disp ' '; 214 fmt='%16.6f %12.6f %8.2f'; 215 results=[beta se tstat]; 216 for i=1:K; 217 fprintf(indv(i,:)); 218 fprintf(1,fmt,results(i,:)); 219 fprintf('\n'); 220 end 221 disp ' '; 222 disp ' '; 223 disp '============================================================================='; 224 end; % if PRINT==0 225 226 K=K+prevest; % Return Correct # of Estimated Pars 227 228 %%%%%%%%%%%%%%% End of OLS %%%%%%%%%%%%%%%% 229 230 231 if FE==0; 232 % Now let's also do a Breusch-Pagan (1980) test of H0: sig2u=0, which 233 % follows Greene (1990), page 491ff.==> Use OLS pooled residuals! 234 235 T=NT/NN; 236 [LM pval]=BPagan(reshape(u,T,NN)); 237 % num=sum((sum(reshape(u,T,NN)).^2)'); 238 % den=u'*u; 239 % LM=NN*T/2/(T-1)*(num/den-1)^2; 240 % pval=1-chicdf(LM,1); 241 % fprintf('B/Pagan Test of H0:Var(ui)=0: LM=%6.2f P-Value=%4.2f\n',... 242 % [LM pval]); 243 244 end; 245 246 247 if FE==4; % Hausman Test, p 495 Greene 248 249 % Now a Hausman test of the orthogonality of random effects 250 if theta>0; 251 Kp=K-prevest; 252 Sigma=vcvW-vcv(2:Kp,2:Kp); 253 bdiff=bW-beta(2:Kp); 254 W=bdiff'*inv(Sigma)*bdiff; 255 pval=1-chicdf(W,Kp-1); 256 fprintf('Hausman Test of E(e|X)=0: W= %6.2f P-Value=%4.2f\n',[W pval]); 257 else; 258 disp '(Theta is Negative -- No Hausman Test conducted)'; 259 end; 260 261 disp ' '; 262 end; 263 NN=N; 264 end;