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