%clean and plot the log wage posterior results
clear;
clc;
randn('seed',sum(100*clock));

load wages.raw;

wage = wages(:,1);
ed  = wages(:,2);
able = wages(:,3);
D = round(.5 + .5*sign(able));

y = log(wage);
nobs = length(y);
x = [ones(nobs,1) ed able ];
k = size(x,2);

%prior parameters
a=3;
b = 2.5;

V_beta = 4*eye(k);
invV_beta = inv(V_beta);

%define some stuff useful for implementing the simulator
points_high = find(D==1);
points_low = find(D==0);
X_high = x(points_high,:);
X_low = x(points_low,:);
y_high = y(points_high);
y_low = y(points_low);
n_high = length(points_high);
n_low = length(points_low);


iter = 1500;
burn = 100;
betas_final = zeros(iter-burn,k);
sighigh_final = zeros(iter-burn,1);
siglow_final = zeros(iter-burn,1);


invSigma = eye(nobs); %full inverse covariance matrix initial condition
for i = 1:iter 
    %sample beta
    D_beta = inv(x'*invSigma*x + invV_beta);
    d_beta = x'*invSigma*y;
    H_beta = chol(D_beta);
    betas = D_beta*d_beta + H_beta'*randn(k,1);
    
    %sample sighigh
    resids_high = y_high - X_high*betas;
    sighigh = invgamrnd( (n_high/2) + a, inv( inv(b) + .5*sum(resids_high.^2)),1,1);
    
    
    %sample siglow
    resids_low = y_low - X_low*betas;
    siglow = invgamrnd( (n_low/2) + a, inv( inv(b) + .5*sum(resids_low.^2)),1,1);
    
    %construct Sigma given both variance parameters
    variance_vec = zeros(nobs,1);
    variance_vec(points_high) = sighigh;
    variance_vec(points_low) = siglow;
    invSigma = diag(1./variance_vec);

    if i > burn;
        betas_final(i-burn,:) = betas';
        sighigh_final(i-burn,1) = sighigh;
        siglow_final(i-burn,1) = siglow;
    end;
end;

disp('post means and std of betas');
[mean(betas_final); std(betas_final)]

disp('posterior mean and std of high ability variance');
[mean(sighigh_final); std(sighigh_final)]

disp('post mean and std of low ability variance')
[mean(siglow_final); std(siglow_final)]

disp('Probability that low ability variance is larger than high ability variance')
mean(.5 + .5*sign(siglow_final - sighigh_final))
timesbigger = (siglow_final)./sighigh_final;
hist(timesbigger,20);