%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ordered Probit a la Poirier & Tobias Book chapter
% i.e. re-parameterized a la Nandram & Chen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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\aaec6564\mlab\worksp\paddlers;
% This is a survey implemented at Lake Tahoe (CA) in 2008. We intercepted
% kayakers and other paddlers at campgrounds and beaches to ask them about
% their paddling experience at Lake Tahoe. The policy background was a
% planned "water trail" (map with suggested landing points around the
% lake).
%
% Among others, we gave them a few Likert-scale questions for the main
% reasons they chose Tahoe (5-point scale from "not at all important" to
% "very important."
%
% A key question was the importance of beach access (many local residents
% were concerned that visiting paddlers would trespass onto their private
% beaches).
%
% The data include 279 observations. Some respondents were visitors, others
% were locals.
% Variable Obs Mean Std. Dev. Min Max
%
% 1 survey_id
% 2 yrs_padd years of paddling experience
% 3 scenic reason for paddling at Tahoe: scenic beauty
% ( 0-4 Likert: not at all important to very important)
% 4 beaches reason for paddling at Tahoe: boat access to beaches
% 5 wildlife reason for paddling at Tahoe: view wildlife
% 6 solitude reason for paddling at Tahoe: experience solitude
% 7 family reason for paddling at Tahoe: be with family / friends
% 8 workout reason for paddling at Tahoe: physical challenge /exercise
% 9 gender 0=male, 1=female
% 10 age age in years
% 11 resident 0 = visitor, 1 = resident of Tahoe Basin
% 12 inc000 income in $1000
% 13 kids 1=kids present in household
% 14 college 1=collge eductaion or higher
y=data(:,4); %let's take a closer look at "beaches"
n=length(y);
f0=find(y==0);
f1=find(y==1);
f2=find(y==2);
f3=find(y==3);
f4=find(y==4);
p0=length(f0)/n;
p1=length(f1)/n;
p2=length(f2)/n;
p3=length(f3)/n;
p4=length(f4)/n;
prop=[p0;p1;p2;p3;p4];
X=[ones(n,1) data(:,2) data(:,9:14)];
k=size(X,2);
% contents of X
% 1 constant
% 2 yrs_padd years of paddling experience
% 3 gender 0=male, 1=female
% 4 age age in years
% 5 resident 0 = visitor, 1 = resident of Tahoe Basin
% 6 inc000 income in $1000
% 7 kids 1=kids present in household
% 8 college 1=collge education or higher
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=15000; % burn-ins
r2=5000; % keepers
R=r1+r2;
nc=3; %number of cutoff points to be estimated
% 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;
% elements for del2 (squared inverse of highest cutoff (c3), but also new
% variance of regression error)
eta0=1/2;
kappa0=1/2; % diffuse prior shape and scale
del2draw=s2; % use OLS variance as starting draw for GS (Tobias uses "1" in his code)
% Priors and Metropolis Hastings elements for cutoffs
Vc=10;
stdc1=0.1;
stdc2=0.1;
c1draw=1/4;
c2draw=1/2; % starting draws for re-parameterized thresholds
[betamat,cmat,acceptmat]=gs_OP(X,y,k,n,nc,r1,r2,mu0,V0,betadraw,eta0,kappa0,del2draw,...
Vc,stdc1,stdc2,c1draw,c2draw,f0,f1,f2,f3,f4);
'GS done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convergence Checks - Geweke's CD diagnostic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open log file
[fid]=fopen('c:\klaus\aaec6564\mlab\logs\mod11_op_paddlers.txt','w');
if fid==-1;
warning('File could not be opened');
else;
disp('File opened successfully');
end;
fprintf(fid,'sample proportions of Likert scores =\t%6.3f \n',prop);
fprintf(fid,'\n');
% put all draws together & run diagnostics
allmat=[betamat;cmat];
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,'acceptance rates in MH part =\t%6.3f \n',acceptmat);
% 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');
%c stuff
out=kdiag(k+1:end,:)';
fprintf(fid,'Output table for thresholds \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\mod11_op_paddlers betamat cmat;
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;