%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PS2 Q1a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;

% 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

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

% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
% load c:\klaus\AAEC6564\mlab\worksp\wetlands.txt;
% data = wetlands;
% clear wetlands;
% save c:\klaus\AAEC6564\mlab\worksp\wetlands;
% 'ok'
% break

load c:\klaus\AAEC6564\mlab\worksp\wetlands;

% variable list

%1  study_id        study ID
%2  obs_id          w/in study observation ID
%3  users           percentage of active wetland users in sample            
%4  lnwtp           log of annual wtp, in 2006 dollars             
%5  lninc           log of mean (if reported for sample) or median 
%                   (if taken from census) HH income in 2006 dollars                          
%6  lnacres         log of wetland acres                    
%7  lnac2           log of wetland acres, squared  

acres=exp(data(:,6));
ac000=acres/1000;

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

X=[ones(n,1) data(:,5) data(:,3) ac000];
% % contents of X:
% 
% 1     constant                
% 2     log(income)
% 3     users
% 4     wetland acres, in units of 1000

k=size(X,2); 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=5000; % 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
V0=eye(k)*100; % 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]=...
    gs_normal_independent(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw);

'GS done'


[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod2_wetlands.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');

% 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');


save c:\klaus\AAEC6564\mlab\worksp\mod2_wetlands betamat sig2mat;

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;








