close all;
clear;
clc;

%Econ 678: Take-Home Exam
%This little pup generates some normal data from a linear process with
%normal errors and then estimates a mis-specified spatial auto-
%regressive model. 

%Initializing the paramaters of the true model
beta = 3;
gamma = 2;

%Initializing the parameters used for simulating the data. 
%Use the following order throughout: [x , z]^T
%The number of data samples 
m =10000;
%The number of observations in each sample
n = 1000;
randn('seed',123456789);
normrnd('seed',123456789);
%Now define the errors u
u = normrnd(0,9,n,m);
%Set the means for x and z
mu = [0;0];
%Change the covar between x and z
covar = [0 0.1 0.2 0.4 0.6 0.8 0.9 0.99];
%Predefine the matrix that will hold the lambda_mle's 
lambda_mle_m = 99*ones(m,size(covar,2));

for k = 1:size(covar,2);
Sigma = [1,covar(1,k);covar(1,k),1];
%mvnrnd generates each row of R using the corresponding row of mu.
%For us this means
xz = mvnrnd(mu,Sigma,n);  
x=xz(:,1);
z=xz(:,2);

%Make the y's.  
xbig = repmat(x,1,m);
zbig = repmat(z,1,m);
ybig = beta*xbig + gamma*zbig + u;
y = ybig(:,1);

%Create W matrix
%First get the differences
Wdif = repmat(z,1,n)-repmat(z',n,1);
%Second square them
Wdif2 = Wdif.^2+eye(n,n);
%Third invert them element by element and put zeros on the diagonal
W = 1./Wdif2-eye(n,n);
%Fourth, normalize the matrix to have row sums of one.
Wnorm = mat2gray(W);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximizing the lambda concentrated log-likelihood 
x0 = 0.5; %initial guess for the fmin function
%Returns the lambda that minimizes the function. 

for i=1:m
lambda_mle_m(i,k) = fminsearch(@(theta) lambdacon(theta,ybig(:,i),x,Wnorm),x0);
end

end %End the loop for different covars
save('MonteCarlo_out.mat');