%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PS1 Q3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Application: home sales in CT

% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
[your turn - uncomment the following 7 lines (right-click, "uncomment"), run them by right-clicking & 
    "Evaluate selection", then comment them again.
% load c:\klaus\AAEC6564\mlab\worksp\CTsales.raw;
% data=CTsales;
% clear CTsales;
% save c:\klaus\AAEC6564\mlab\worksp\CTsales data;
% 
% 'ok'
% return

load c:\klaus\AAEC6564\mlab\worksp\CTsales;
%7471 observations on residential home sales along the Connecticut coast
%between 2002 and 2014; avg. price = 476,000; max = 1.4 million.

% 1     DQid                            original SR_PROPERTY_ID
% 2     saleyr                           year of sale
% 3     salemonth                        month of sale
% 4     saleday                          calendar day of sale
% 5     saledateE                        sale date in elapsed days since 1/1/1960
% 6     price000                         sale price in 1000's of 2014 dollars
% 7     wf                               coastal waterfront property, 0.05 mile rule
% 8     hardshore                        1= nearest shoreline is hardened (= seawall or equivalent)
% 9     age                              age of structure, years
% 10    sqft00                           square footage, in 100's
% 11    lot000                           lot size, in 1000 sqft
% 12    bedrooms                        
% 13    bathrooms                        total number of baths
% 14    elev10                           elevation in meters at 10 meter resolution
% 15    ISMi                             distance to nearest Interstate, miles
% 16    PAMi                             distance to nearest principal artery, miles
% 17    beaMi                            distance to nearest beach, miles
% 18    hidMi                            distance to nearest high-density development, miles
% 19    coaestMi                         distance to nearest coast or estuary, miles
% 20    reslkMi                          miles to nearest lake, pond, or reservoir
% 21    rivMi                            Distance in miles to nearest stream or river
% 22    ag10                             acreage of ag land w/in 1000m, most current
% 23    ind10                            acreage of ind. land w/in 1000m, most current
% 24    op10                             acreage of open land w/in 1000m, most current
% 25    SY12                             saleyr== 2002.0000
% 26    SY13                             saleyr== 2003.0000
% 27    SY14                             saleyr== 2004.0000
% 28    SY15                             saleyr== 2005.0000
% 29    SY16                             saleyr== 2006.0000
% 30    SY17                             saleyr== 2007.0000
% 31    SY18                             saleyr== 2008.0000
% 32    SY19                             saleyr== 2009.0000
% 33    SY20                             saleyr== 2010.0000
% 34    SY21                             saleyr== 2011.0000
% 35    SY22                             saleyr== 2012.0000
% 36    SY23                             saleyr== 2013.0000
% 37    SY24                             saleyr== 2014.0000

y=data(:,6);
n=length(y);

X=[ones(n,1) data(:,7:24) data(:,26:end)];
k=size(X,2);

% Contents of X:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1     constant term
% 2     wf                               coastal waterfront property, 0.05 mile rule
% 3     hardshore                        1= nearest shoreline is hardened (= seawall or equivalent)
% 4     age                              age of structure, years
% 5     sqft00                           square footage, in 100's
% 6     lot000                           lot size, in 1000 sqft
% 7     bedrooms                        
% 8     bathrooms                        total number of baths
% 9     elev10                           elevation in meters at 10 meter resolution
% 10    ISMi                             distance to nearest Interstate, miles
% 11    PAMi                             distance to nearest principal artery, miles
% 12    beaMi                            distance to nearest beach, miles
% 13    hidMi                            distance to nearest high-density development, miles
% 14    coaestMi                         distance to nearest coast or estuary, miles
% 15    reslkMi                          miles to nearest lake, pond, or reservoir
% 16    rivMi                            Distance in miles to nearest stream or river
% 17    ag10                             acreage of ag land w/in 1000m, most current
% 18    ind10                            acreage of ind. land w/in 1000m, most current
% 19    op10                             acreage of open land w/in 1000m, most current
% 20    SY13                             saleyr== 2003.0000
% 21    SY14                             saleyr== 2004.0000
% 22    SY15                             saleyr== 2005.0000
% 23    SY16                             saleyr== 2006.0000
% 24    SY17                             saleyr== 2007.0000
% 25    SY18                             saleyr== 2008.0000
% 26    SY19                             saleyr== 2009.0000
% 27    SY20                             saleyr== 2010.0000
% 28    SY21                             saleyr== 2011.0000
% 29    SY22                             saleyr== 2012.0000
% 30    SY23                             saleyr== 2013.0000
% 31    SY24                             saleyr== 2014.0000


% Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=10000; % burn-in
r2=10000; % keepers
R=r1+r2;

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(n-k);

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
c=100;
V0=c*eye(k); % diffuse prior for varcov of beta
betadraw=bols; %use OLS values as starting draws for Gibbs Sampler

% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
v0=1/2;
tau0=1/2; % diffuse prior shape and scale

[betamat,sig2mat]= ...
    gsPS1Q2(X,y,n,k,r1,r2,mu0,V0,v0,tau0,sig2draw);


'GS done'

[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\ps1q3klaus.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% put all draws together & run diagnostics
allmat=[betamat;sig2mat];
kdiag=klausdiagnostics_greater0(allmat);

fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');

fprintf(fid,'OLS betas \t%6.3f\n',bols); 
fprintf(fid,'OLS s2 \t%6.3f\n',s2); 
fprintf(fid,'\n');
fprintf(fid,'\n');

% beta stuff
out=kdiag(1:k,:)';
fprintf(fid,'Output table for betas \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

%sig2 stuff
out=kdiag(k+1,:)';
fprintf(fid,'Output table for sig2 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


%focus on selected marginal effects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wf=betamat(2,:)';
hardshore=betamat(3,:)';
bedrooms= [your turn]
ISmiles= [your turn]

[f1,x1]=ksdensity(wf,'kernel','epanechnikov','npoints',500);
[f2,x2]=ksdensity(hardshore,'kernel','epanechnikov','npoints',500);
[f3,x3]=ksdensity(bedrooms,'kernel','epanechnikov','npoints',500);
[f4,x4]=ksdensity(ISmiles,'kernel','epanechnikov','npoints',500);

figure(1);

subplot(2,2,1);
plot(x1,f1,'-k');
title('posterior for waterfront effect');
set(gca,'Xlim',[-10 150]);
set(gca,'Ylim',[0 0.05]);
xline(0,'--r'); %add red dashed reference line at zero

subplot(2,2,2);
plot(x2,f2,'-k');
title('posterior for hardened shoreline');
set(gca,'Xlim',[your turn - find limits that look OK and include zero]);
set(gca,'Ylim',[yourturn]);
xline(0,'--r'); %add red dashed reference line at zero

subplot(2,2,3);
plot(x3,f3,'-k');
title('posterior for bedrooms');
[your turn]

subplot(2,2,4);
[your turn]

save c:\klaus\AAEC6564\mlab\worksp\ps3klaus betamat sig2mat X y;

finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);


st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;










