% Program to simulate Nakamura and Nakamura (N&N), 1985, Table 1
% The main loop simulates the regressions in N&N for different sample
% sizes (here, 20 100 250). The number of monte carlos is controlled by
% nsim (here 200,000).
% At the end of the main loop is a short block of code that does monte
% carlos by simulating Imhof's method over different deltas.
%
% Copyright, PCR 2016. I give no warranty for the code.
clear; close all; clc;
% RNG Initialization 
state = 16483;         % Yes, prior RNG; but allows duplication.
rand('state',state);
randn('state',state);
% N     = Number of endogenous variables
% K     = Number of exogenous variables
  N=2; K=6;
% Nsim  = Number of monte carlo simulations 
  nsim=200000;
% Table_res = Results matrix
% Tvals = Sample Sizes
Table_res=[];
% R_2   = Population R^2 values for reduced form
% rho2  = Correlation coefficient
% B     = See text

int_result4=[];
% Main loop does different sample sizes
for Tvals=[20 100 250];
    % This loop does different population R^2s
    for R_2=[.2 .5 .8];
        % This loop does different error correlations 
        for rho2=[0 .1 .5];
            % This loop does different B parameters (see N&N 1985)
            for B=[ .2 .8];
                % Variance parameters. Assumes x's are standard normal. See text
                lambd(2)=sqrt(5*rho2*(1-R_2)/R_2);
                lambd(3)=sqrt(5*(1-rho2)*(1-R_2)/R_2);
                lambd(1)=10;
            if lambd(2) > 0 
                lambd(1)=B*(3+lambd(2)^2+lambd(3)^2)/lambd(2);
            end
     T=Tvals;
     y_1=[];
     Y=[];
% omeg_x = Variance-Covariance matrix of exogenous variables
     omeg_x=eye(K-1,K-1);  
     omeg_x_half=chol(omeg_x);    % Cholesky Decomposition
% x = stochastic random variables
% X = [1 x] includes a constant term
     x=normrnd(zeros(T,K-1),1);
     x=x*omeg_x_half;
     X2=x;  
% omeg_e = Variance-Covariance matrix of structural eq errors
     omeg_e=eye(N,N);
     omeg_e_half=chol(omeg_e);
     omeg_e(2,2)=(lambd(2)^2+lambd(3)^2);
     omeg_e(1,1)=lambd(1)^2;
     % gam is delta ib=n paper Cov(eta1, V2)/Var(V2)
     gam=lambd(1)*lambd(2)/(lambd(2)^2+lambd(3)^2);
     % Structural parameters
     bvec=[1 1 1 0 ]'; pi1=[1 1 0 0 0]';pi2=[1 1 1 1 1]';
     bvec2=[1 1 1 gam']';
     % F test restriction matrix/vector
     R=[0 0 0 1]; r=[0];
     res=[];
     % Simulations of N&N (1985) regressions 
     for j=1:nsim;
        x=normrnd(zeros(T,K-1),1); X=x;
        % The errors
        e=normrnd(zeros(T,N),1);
        e=[lambd(1)*e(:,1) (lambd(2)*e(:,1)+lambd(3)*e(:,2))];
        % Dependent variables
        Y(:,2)=X*pi2+e(:,2);
        Y(:,1)=Y(:,2)+X*pi1+e(:,1);
        % Regressors
        Vhat=Y(:,2)-X*inv(X'*X)*X'*Y(:,2);
        Z=[X(:,1:2) Y(:,2) Vhat];
        y_1=Z*bvec+e(:,1);
        Z1=Z(:,1:end-1);    % No Vhat
        W=Z*inv(Z'*Z)*R';
        bhat=inv(Z'*Z)*Z'*y_1;
        btilde=bhat+inv(Z'*Z)*R'*inv(R*inv(Z'*Z)*R')*(r-R*bhat);
        btilde2=inv(Z1'*Z1)*Z1'*y_1;
        num=(r-R*bhat)'*inv(R*inv(Z'*Z)*R')*(r-R*bhat);
        num2=(y_1-Z*btilde)'*(y_1-Z*btilde)-(y_1-Z*bhat)'*(y_1-Z*bhat);
        denom=(y_1-Z*bhat)'*(y_1-Z*bhat);
        mu=(e(:,2)-Vhat)*gam;
        delt1=(R*bvec2-r+W'*mu)'*inv(W'*W)*(R*bvec2-r+W'*mu)/(omeg_e(1,1)-gam'*omeg_e(2:N,2:N)*gam);
        delt2=(mu'*mu-mu'*Z*inv(Z'*Z)*Z'*mu)/...
        (omeg_e(1,1)-gam'*omeg_e(2:N,2:N)*gam);
        % res stores rsq plus Deltas eqn (33)
        res(j,:)=[num/(num+denom) delt1 delt2];
   % Prints out progress every 50K iterations
   if mod(j,50000)==0
       [j Tvals R_2 rho2 B] 
   end
end;

rsq=sortrows(res(:,1),1);
f_crit=finv(.95,1,T-4);
rsq_crit=(f_crit/(T-4))/(1+(f_crit/(T-4)));
indic=rsq > rsq_crit;


% My Code for simulating the imhof integral 
delt=[0 0]';
multip=[size(R,1) T-size(bvec,1)];
int_result=[];
u_bnd=30;  % Upper bound of integral (set to something large)
arg=0;
kk=0;
int_ans=[];int_ans2=[];
for reps=[10 50 1000];    % Number of Imhof integrals averaged (Averaging over deltas in delt)
    int_result3=zeros(1,reps);
    lambdv=[1-rsq_crit; -rsq_crit];
    for k=1:reps
      % Note res contains r2 delt1 delt2 for each simulated data set
      delt=[res(k,2:3)]';
      int_result3(1,k)=.5-(1/pi)*quad(@imhof_integrand,1e-6,u_bnd,1e-6,0,arg,lambdv,multip,delt);
    end;
  % Stores average and SE of average of Imhof integrals
  int_ans=[int_ans mean(1-int_result3')];
  int_ans2=[int_ans2 std(1-int_result3')/sqrt(reps)];
end;
% Stores results: First three cols are N&N 1985 cols; 4) N; 5) F Crti;
%   6) R^2 crit; 7)
Table_res=[Table_res; [R_2 B rho2 T f_crit rsq_crit sum(indic)/nsim int_ans int_ans2]]

end;
end;
end;

end; %tvals loop
sprintf('%3.1f %3.1f %3.1f %4.0f %7.5f %7.5f %7.5f %7.5f\n  %7.5f %7.5f %7.5f\n',[Table_res(:,1:4) Table_res(:,7:13)]')
