% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Poisson model using TCW fishing data
% In reality, each angler reported trips to 12 sites, which would suggest a
% hierarchical model
% However, for this simple example, we will only use site #4 (vicinity of Reno)
% and treat each observation as independent
% n=551
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;  % start stop watch

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

% load c:\klaus\mudsnail2\mlab\worksp\fromSTATA.raw;
% data=fromSTATA;
% clear fromSTATA;
% save c:\klaus\mudsnail2\mlab\worksp\TCWfishing;
% 'ok'
% break

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

% Contents of data (sorted by individuals and sites)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1 id         person id  (1-551)                    
%2 site       site id   (1-12)                
%3 trips      2004 trips to each site                    
%4 gender     1=female
%5 age                      
%6 fly        1=flyfishing only                      
%7 spin       1=spin casting only (baseline = both)                    
%8 edu        years of formal education
%9 lninc      log of annual HH income
%10 fshyrs    years of fishing
%11 kids      1=kids in HH
%12 tc        travel cost (see notes / STATA do file)                
%13 cr        1= catch and realease at site
%14 quant     1= quantity restrictions (other than cr) at site
%15 special   1= special regulations at site (comprises cr + quant)
%16 winter    1= winter closure (last Sat. Apr. - Nov.15)
%recode gender 2=0
f=find(data(:,4)==2);
data(f,4)=0;

% keep only Truckee #4 - essentially Reno/Sparks area
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=find(data(:,2)==4);
data=data(f,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y=data(:,3);
N=length(y); 

price=data(:,12); %travel cost to access site
inc000=(exp(data(:,9)))/1000;


X=[ones(n,1) price data(:,4:5) data(:,5).^2 data(:,6:7) inc000 ...
   data(:,10:11)];
k=size(X,2);

%Contents of X
%%%%%%%%%%%%%%%%%%%%%%%%%%%               
%1  constant
%2  tc         (dollars)
%14 gender     1=female
%15 age
%16 age2
%17 fly        1=flyfishing only                      
%18 spin       1=spin casting only (baseline = both)
%19 inc000     log of annual HH income
%20 fshyrs     years of fishing
%21 kids       1=kids in HH  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=10000; % burn-ins - let the CD diagnostic be your guide...
r2=10000; % keepers
R=r1+r2;

% generic OLS
bols=inv(X'*X)*X'*y;

% elements for beta (fixed coefficients)
mub0=zeros(k,1); 
Vb0=eye(k)*100; 

tau=1; %tuner for variance in t-distribution
v=10; %degrees of freedom in t-distribution
betadraw = bols;

% 'ok'
% return

[betamat,acceptmat]=gs_Poisson_ic(y,X,k,N,mub0,Vb0,tau,v,betadraw,r1,r2);

'GS done'


% put all draws together
allmat=[betamat];%
kdiag=klausdiagnostics_greater0(allmat);

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

fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'tuner for mvt variance in MH part=\t%6.3f \n',tau);
fprintf(fid,'DoF in mvt=\t%6.3f \n',v);
fprintf(fid,'acceptance rates in MH part =\t%6.3f \n',acceptmat);
fprintf(fid,'\n');

% beta stuff
out=kdiag'; 
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');

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;

save c:\klaus\AAEC6564\mlab\worksp\mod7_Poisson_ic betamat;

acceptmat