%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare convergence plots under different blocking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Re-run Efficient 2-block model, keep burn-ins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load c:\klaus\AAEC6564\mlab\worksp\labordata;
%We will use Mroz's (1987) labor data described in Greene (5th edition),
%p. 51
% Contents of Data (columns)
%1 LFP = A dummy variable = 1 if woman worked in 1975, else 0
%2 WHRS = Wife's hours of work in 1975
%3 KL6 = Number of children less than 6 years old in household
%4 K618 = Number of children between ages 6 and 18 in household
%5 WA = Wife's age
%6 WE = Wife's educational attainment, in years
%7 WW = Wife's average hourly earnings, in 1975 dollars
%8 RPWG = Wife's wage reported at the time of the 1976 interview (not = 1975 estimated wage)
%9 HHRS = Husband's hours worked in 1975
%10 HA = Husband's age
%11 HE = Husband's educational attainment, in years
%12 HW = Husband's wage, in 1975 dollars
%13 FAMINC = Family income, in 1975 dollars
%14 WMED = Wife's mother's educational attainment, in years
%15 WFED = Wife's father's educational attainment, in years
%16 UN = Unemployment rate in county of residence, in percentage points.
%17 CIT = Dummy variable = 1 if live in large city (SMSA), else 0
%18 AX = Actual years of wife's previous labor market experience
f=find(data(:,1)==0); % we only include women that worked in 1975
data(f,:)=[];
y=log(data(:,2).*data(:,7)); %dependent variable = log(hrs worked* hourly wage)
KL6=data(:,3); %number of children under 6
K618=data(:,4); % Number of children between ages 6 and 18 in household
age=data(:,5);
age2=age.^2;
edu=data(:,6); %years of education
n=length(y);
X=[ones(n,1) age age2 edu KL6 K618];
k=size(X,2);
r1=5000; % burn-in
r2=10000; % keepers
R=r1+r2;
% 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; %use OLS values as starting draws for Gibbs Sampler
% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
v0=1/2;
tau0=1/2; % diffuse prior shape and scale
[betamat,sig2mat]=...
gs_normal_independent_keepall(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw);
betamat1=betamat;
sig2mat1=sig2mat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Re-run inefficient 4-block model, keep burn-ins
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X1=X(:,1:2);
X2=X(:,3:4);
X3=X(:,5:6);
k1=size(X1,2);
k2=size(X2,2);
k3=size(X3,2);
k=k1+k2+k3;
% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(n-k);
% elements for beta1
mu01=zeros(k1,1);
V01=eye(k1)*100;
beta1draw=bols(1:k1);
% elements for beta2
mu02=zeros(k2,1);
V02=eye(k2)*100;
beta2draw=bols(k1+1:k1+k2);
% elements for beta3
mu03=zeros(k3,1);
V03=eye(k3)*100;
beta3draw=bols(k1+k2+1:end);
% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
v0=2.1;
tau0=20; % diffuse prior shape and scale
[beta1mat,beta2mat,beta3mat,sig2mat]=...
gs_normal_blocked_keepall(X,X1,X2,X3,y,n,k,k1,k2,k3,r1,r2,mu01,V01,beta1draw,...
mu02,V02,beta2draw,mu03,V03,beta3draw,v0,tau0,sig2draw);
betamat2=[beta1mat;beta2mat;beta3mat];
sig2mat2=sig2mat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare convergence for "age"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=(1:1:R);
figure(1);
subplot(2,1,1);
age=betamat1(2,:);
plot(x,age,'b');
axis([min(x) max(x) min(age) max(age)]);
title('age, 2-block sampler');
xlabel('iteration');
ylabel('drawn values')
subplot(2,1,2);
age=betamat2(2,:);
plot(x,age,'r');
axis([min(x) max(x) min(age) max(age)]);
title('age, 4-block sampler');
xlabel('iteration');
ylabel('drawn values')
figure(2);
subplot(2,1,1);
edu=betamat1(4,:);
plot(x,edu,'b');
axis([min(x) max(x) min(edu) max(edu)]);
title('education, 2-block sampler');
xlabel('iteration');
ylabel('drawn values')
subplot(2,1,2);
edu=betamat2(4,:);
plot(x,edu,'r');
axis([min(x) max(x) min(edu) max(edu)]);
title('education, 4-block sampler');
xlabel('iteration');
ylabel('drawn values')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Show bivariate convergence plot for age, age2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
subplot(2,1,1);
age=betamat1(2,:);
age2=betamat1(3,:);
plot(age,age2,'b');
title('age2 against age, efficient sampler');
xlabel('age');
ylabel('age2')
set(gca,'Xlim',[-0.2 0.5]);
set(gca,'Ylim',[-0.01 0.01]);
subplot(2,1,2);
age=betamat2(2,:);
age2=betamat2(3,:);
plot(age,age2,'r');
title('age2 (block 2) against age (block 1), inefficient sampler');
xlabel('age');
ylabel('age2')
set(gca,'Xlim',[-0.2 0.5]);
set(gca,'Ylim',[-0.01 0.01]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Show bivariate convergence plot for edu, K618
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);
subplot(2,1,1);
edu=betamat1(4,:);
K618=betamat1(6,:);
plot(edu,K618,'b');
title('edu against K618, efficient sampler');
xlabel('edu');
ylabel('K618')
set(gca,'Xlim',[-0.05 0.2]);
set(gca,'Ylim',[-0.4 0]);
subplot(2,1,2);
edu=betamat2(4,:);
K618=betamat2(6,:);
plot(edu,K618,'r');
title('edu (block 2) against K618 (block 3), inefficient sampler');
xlabel('edu');
ylabel('K618')
set(gca,'Xlim',[-0.05 0.2]);
set(gca,'Ylim',[-0.4 0]);