% suggested solutions to hw10, question 3

sigmau = 1;
sigmav = 1;
sigmaz = 1;
beta = 1;
lambda = 1;
n = 100;

J = 1000;
gamma = 2;

biv_draws_a = NaN*ones(J,1);
for j=1:J,
	u = sigmau*randn(n,1);
        v = sigmav*randn(n,1);
        z = sigmaz*randn(n,1);
        x = gamma*z + gamma*v + lambda*u;
        y = beta*x + u;
        
        biv_draws_a(j,1) = (z'*y)/(z'*x);	
end;
% remark: since the main monte carlo routine is      
% repeated, we could have written it as a function 
% and saved some typing.

hist(biv_draws_a,30)
title('\gamma = 2')
mean(biv_draws_a)
var(biv_draws_a)

gamma = .5;
biv_draws_b = NaN*ones(J,1);
for j=1:J,
	u = sigmau*randn(n,1);
        v = sigmav*randn(n,1);
        z = sigmaz*randn(n,1);
        x = gamma*z + gamma*v + lambda*u;
        y = beta*x + u;
        
        biv_draws_b(j,1) = (z'*y)/(z'*x);	
end;

truncation = (biv_draws_b>=-8)&(biv_draws_b<=8);
% remark: 'truncation' is a logical vector, which is used
%          below to select only observations in [-8,8]

figure
hist(biv_draws_b(truncation),30 )
title('\gamma = .5')
mean(biv_draws_b)
var(biv_draws_b)


gamma = .1;
biv_draws_c = NaN*ones(J,1);
for j=1:J,
	u = sigmau*randn(n,1);
        v = sigmav*randn(n,1);
        z = sigmaz*randn(n,1);
        x = gamma*z + gamma*v + lambda*u;
        y = beta*x + u;
        
        biv_draws_c(j,1) = (z'*y)/(z'*x);	
end;

truncation = (biv_draws_c>=-8)&(biv_draws_c<=8);

figure
hist(biv_draws_c(truncation),30)
title('\gamma = .1')
mean(biv_draws_c)
var(biv_draws_c)

% end of hw10_3.m