% PRETEST ESTIMATORS clc; % Define Parameters and Independent Variables n = 50; loops = 2000; x0 = ones(n,1); x2 = randn(n,1); x1 = 10*rand(n,1)-x2; xmat = [x0,x1,x2]; invxx = inv(xmat'*xmat); b = [0.5; 1; 2]; k = size(b,1); %Simulate Artificial Data y = zeros(n,loops); for i = (1:1:loops); sig = 10; y(:,i) = xmat*b + sig*randn(n,1); end; % Unrestricted Estimator urb = zeros(k,loops); err = zeros(n,loops); s2 = zeros(loops,1); for i = (1:1:loops) urb(:,i) = invxx*(xmat'*y(:,i)); err(:,i) = y(:,i) - xmat*urb(:,i); s2(i,1) = (err(:,i)'*err(:,i))/(n - k); end % Restricted Estimator rb = zeros(k-1,loops); rxmat = xmat(:,1:2); for i = (1:1:loops) rb(:,i) = inv(rxmat'*rxmat)*(rxmat'*y(:,i)); end % Pretest Estimator tcrit = 2.0; seurb = zeros(loops,1); tstat = zeros(loops,1); ptb = zeros(loops,1); for i = (1:1:loops) seurb(i,1) = s2(i,1)*invxx(3,3); tstat(i,1) = urb(k,i)/sqrt(seurb(i,1)); if abs(tstat(i,1)) < tcrit; ptb(i,1) = rb(k-1,i); else; ptb(i,1) = urb(k-1,i); end end % Mean Squared Error Criterion urmse = zeros(loops,1); rmse = zeros(loops,1); ptmse = zeros(loops,1); for i = (1:1:loops) urmse(i,1) = (urb(2,i) - b(2,1)).^2; rmse(i,1) = (rb(2,i) - b(2,1)).^2; ptmse(i,1) = (ptb(i,1) - b(2,1)).^2; end avgurmse = mean(urmse); avgrmse = mean(rmse); avgptmse = mean(ptmse); disp(['Unrestricted Average MSE = ' num2str(avgurmse')]); disp(' '); disp(['Restricted Average MSE = ' num2str(avgrmse')]); disp(' '); disp(['Pretest Average MSE = ' num2str(avgptmse')]); subplot(3,1,1); hist(urb(2,:),50); title('Unrestricted'); set(gca,'XLim',[-3,5]); subplot(3,1,2); hist(rb(2,:),50); title('Restricted'); set(gca,'XLim',[-3,5]); subplot(3,1,3); hist(ptb(:,1),50); title('Pretest'); set(gca,'XLim',[-3,5]);