% PROGRAM NAME: houseprices.m
% From problem 14.53 on p.475 of Miller & Miller
% DATA:
x = [70.3 102 62.5 74.8 57.9 81.6 110.4 88];
y = [114.4 169.3 106.2 125 99.8 132.1 174.2 143.5];
XY = [x' y'];
% Yhat = mX + b model parameter estimates
CXY = cov(XY);
m = CXY(1,2)/CXY(1,1);
MXY = mean(XY);
b = MXY(2) - m*MXY(1);
yhat = m*x + b; yhat = yhat';
% Correlation coefficient & coefficient of determination
r = m*(CXY(1,1)/CXY(2,2))^.5;
R2 = r^2;
figure(1)
plot(x,y,'*',x,yhat,'r')
xlabel('Assessed Value ($1000)')
ylabel('Sell Price ($1000)')
title('Scatter Plot of Assessed vs Sell Prices')
grid
pause
%===============================================
% Use simulations to otain pdf of R2
mx = MXY(1); my = MXY(2);
sigx = CXY(1,1)^.5; sigy = CXY(2,2)^.5; covxy = CXY(1,2);
sigw = (sigy^2 - (covxy/sigx)^2 )^.5;
n = 8; nsim = 5000;
X = mx + sigx*randn(n,nsim);
W = sigw*randn(n,nsim);
Y = m*X + b + W;
sX = std(X); sY = std(Y); mX = mean(X); mY = mean(Y);
cXY = zeros(1,nsim);
for k = 1:nsim
    cc = cov(X(:,k),Y(:,k)); 
    cXY(k) = cc(1,2);
end
mvec = cXY./sX.^2;
rvec = mvec.*(sX./sY);
R2vec = rvec.^2;
% Obtain histogram-based pdf of R2hat
db = 0.005; bvec = .9+db/2 : db : 1-db/2;
h = hist(R2vec,bvec);
fR2 = (nsim*db)^-1 * h;
figure(2)
bar(bvec,fR2)
title('Histogram-based pdf Estimate for R2 Based on 5000 Sims')
grid