% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Hierarchical Poisson model using TCW fishing data
% LOG1 version uses LaFrance & Hanemann's (1984) LOG1
% IDS, which is exponetially linear in income
%
% This version runs in parallel using a local cluster & "parfor" loops
% instead of "for" loops in the gibbs sampler
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
% Open log file
[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod7_HP_mudsnail_parallel.txt','w');
if fid==-1;
warning('File could not be opened');
else;
disp('File opened successfully');
end;
N=551;
J=12;
n=J*N;
% 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;
y=data(:,3);
%create site dummies
Dmat=repmat(eye(J),N,1);
Pmat=Dmat.*repmat(data(:,12),1,J);
inc000=(exp(data(:,9)))/1000;
% Truckee=zeros(n,1);
% f=find(data(:,2)<5);
% Truckee(f)=1;
% Carson=zeros(n,1);
% f=find(data(:,2)>4 & data(:,2)<9);
% Carson(f)=1;
% Walker=zeros(n,1);
% f=find(data(:,2)>8);
% Walker(f)=1;
%Contents of X
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1 constant
%2 tc site 1
%3 tc site 2
%4 tc site 3
%5 tc site 4
%6 tc site 5
%7 tc site 6
%8 tc site 7
%9 tc site 8
%10 tc site 9
%11 tc site 10
%12 tc site 11
%13 tc site 12
%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
%22 special 1= special regulations at site (comprises cr + quant)
%23 winter 1= winter closure (last Sat. Apr. - Nov.15)
X=[ones(n,1) Pmat data(:,4:5) data(:,5).^2 data(:,6:7) inc000 ...
data(:,10:11) data(:,15:16)];
k=size(X,2);
k1=k-2;%number of fixed coeffs
k2=2;%number of random coeffs
% create panel cells for each person
Xfcell=mat2cell(X(:,1:k1),J*ones(N,1),k1);
Xrcell=mat2cell(X(:,k1+1:k),J*ones(N,1),k-k1);
ycell=mat2cell(y,J*ones(N,1),1);
Xf=cell2mat(Xfcell);
Xr=cell2mat(Xrcell);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=10000; % burn-ins
r2=10000; % keepers
R=r1+r2;
% generic OLS
bols=inv(X'*X)*X'*y;
% elements for bet (fixed coefficients)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mub0=zeros(k1,1);
Vb0=eye(k1)*10;
taub=1.5; %tuner for variance in t-distribution
vb=8; %degrees of freedom in t-distribution
betadraw = bols(1:k1);
MHalgo.beta='tailored mvt';
% elements for gam (random coefficients)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mug0=zeros(k2,1);
Vg0=eye(k2)*10;
taug=2; %tuner for variance in t-distribution
vg=8; %degrees of freedom in t-distribution
MHalgo.gi='tailored mvt';
%elements for VCOV of RCs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v0=k2; % number of RCs
S0=eye(k2); % So these are the diffusest priors possible
Edraw=iwishrnd(S0,v0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parpool('local',12); %start parallel
% my PC has 12 cores, so I'm using them all.
% You might need to change this to 4,6,8 etc
% depending on your computer
% you can find the number of cores by typing "feature('numCores')" in the
% command window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[betamat,gmat,gimat,Emat,acceptmat]=gs_HPpar(y,ycell,X,Xf,Xr,Xfcell,Xrcell,k1,k2,k,N,J,n,...
mub0,Vb0,taub,vb,betadraw,...
mug0,Vg0,taug,vg,v0,S0,Edraw,...
r1,r2);
'GS done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete(gcp); %terminates parallel pool
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Capture key results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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, betaf=\t%6.3f \n',taub);
fprintf(fid,'DoF in mvt=\t%6.3f \n',vb);
fprintf(fid,'tuner for mvt variance in MH part, betar=\t%6.3f \n',taub);
fprintf(fid,'DoF in mvt=\t%6.3f \n',vb);
fprintf(fid,'acceptance rates in MH part for betaf =\t%6.3f \n',acceptmat(1));
fprintf(fid,'acceptance rates in MH part for betar =\t%6.3f \n',acceptmat(2));
fprintf(fid,'\n');
allmat=[betamat;gmat;Emat];
kdiag=klausdiagnostics_greater0(allmat);
% beta stuff
out=kdiag(1:k1,:)';
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');
% gamma stuff
out=kdiag(k1+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');
%E stuff
out=kdiag(k+1:end,:)';
fprintf(fid,'Output table for E \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_HP_mudsnail_parallel betamat gmat gimat Emat;