% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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=250000; % burn-ins
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;
% %variance elements for candidate-generating density (CGD)
% load c:\klaus\AAEC6564\mlab\worksp\mod7_Poisson_rwc_prep betamat;
c=0.00000001;
Vc=eye(k);
betadraw = bols;
[betamat,acceptmat]=gs_Poisson_rwc(y,X,k,N,mub0,Vb0,c,Vc,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_rwc_prep.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 variance in MH part=\t%6.3f \n',c);
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_fishing_rwc_prep betamat;