%% lsmonte.m
%% A Matlab program to do a simple least squares monte carlo

%%	we can do this with a large or small sample
%%	with normal or not normal errors

%% we can also do this with an accurate or less accurate number of monte carlo trials

rand('seed',0);
randn('seed',0);

n = 100;								% try other examples such as n=20 or n=5000.                               
           							% You can find big differences

trials = 2000;  					% try same thing with 20 trials to show effect of mc
betahat = zeros(trials,1);   % set aside space

x = 16.2 + 2*rand(n,1);			% create "fixed" x
sumxsqr = sum(x.^2)

for i=1:trials
   u = rand(n,1)-0.5; s2=1/12;			% uniform, variance 1/12
                                       % u = randn(n,1); s2=1.0; <-
                                       % try standard normal
   y = 1.0*x + u;
   betahat(i) = x\y;				% peculiar Matlab notation for a regression
end

disp(['Theory says mean is 1.0 and variance is ',num2str(s2/sumxsqr)]);
mean(betahat)
cov(betahat)
hist(betahat,50);
