%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The following estimates the coefficents of a simulataneous equations 
% model using OLS.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

% Set the true parameters and placeholders for results.

T = 1000;                                                                   % Number of observations.

beta0 = 0;                                                                  % Intercept in equation 1.
beta1 = 0.5;                                                                % Coefficient on Y2 in equation 1.
beta2 = 0.3;                                                                % Coefficient on X1 in equation 1.
beta4 = 0;                                                                  % Intercept in equation 2.
beta5 = 0.3;                                                                % Coefficient on Y1 in equation 2.
beta6 = 0.2;                                                                % Coefficient on X1 in equation 2.
beta7 = 1;                                                                  % Coefficient on X2 in equation 2.

sigma1 = 1;                                                                 % Standard deviation of u1.
sigma2 = 2;                                                                 % Standard deviation of u2.

sim = 1000;                                                                 % Number of simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The simultaneous equations model is
%           Y1 = beta0 + beta1*Y2 + beta2*x1 + u1
%           Y2 = beta4 + beta5*Y1 + beta6*x1 + beta7*x2 + u2
%
% This code estimates the coefficients in the first equation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the G^{-1} and B matrices.

G_inv = zeros(2,2);                                                         % G^{-1} matrix.

det_g = 1 - beta1*beta5;                                                    % Determinant of G.

G_inv(1,1) = 1/det_g;
G_inv(1,2) = beta5/det_g;
G_inv(2,1) = beta1/det_g;
G_inv(2,2) = 1/det_g;

B = zeros(3,2);                                                             % B matrix.

B(1,1) = beta0;
B(2,1) = beta2;
B(3,1) = 0;
B(1,2) = beta4;
B(2,2) = beta6;
B(3,2) = beta7;

% Generate the data.

x1 = 2 + ((1:1:T)'/T).*normrnd(0,1,T,1);                                    % X1.
x2 = 3 + 0.5*x1 + sigma1.*normrnd(0,1,T,1);                                 % X2.

X = [ones(T,1) x1 x2];                                                      % X.

B_2sls = zeros(sim, size(X,2));                                             % 2SLS results.
B_ols = zeros(sim, size(X,2));                                              % OLS results.

for s = 1:sim
    
    u1 = normrnd(0,sigma1,T,1);                                             % Residuls for equation 1.
    u2 = normrnd(0,sigma2,T,1) + 0.2*u1;                                    % Residuals for equation 2.    
    
    u = [u1 u2];                                                            % Residuals.

    Y = X*B*G_inv + u;                                                      % System of equations.
    Y1 = Y(:,1);                                                            % Y1.
    Y2 = Y(:,2);                                                            % Y2.
    
    % Two-stage Least Squares: Stage 1
    
    %{ MISSING CODE %}

    % Two-stage Least Squares: Stage 2

    X1_hat = [ones(T,1) Y2_hat x1];       
    B_2sls(s,:) = inv(X1_hat'*X1_hat)*X1_hat'*Y1;                           % 2SLS estimates.

    % OLS esimation
    
    X1_ols = [ones(T,1) Y2 x1];    
    B_ols(s,:) = inv(X1_ols'*X1_ols)*X1_ols'*Y1;                            % OLS estimates.

end

% Display the results of the last estimation.

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: 2SLS estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')
    
    %{ MISSING CODE %}

fprintf('--------------------------------------------------------------\n')
fprintf('Empirical Results: OLS estimates for Equation 1\n')
fprintf('--------------------------------------------------------------\n\n')
    
    %{ MISSING CODE %}
    
    