Code covered by the BSD License  

Highlights from
Toolkit on Econometrics and Economics Teaching

image thumbnail

Toolkit on Econometrics and Economics Teaching

by

 

Many MATLAB routines related to econometrics, statistics and introductory economics teaching.

AGGREGATE_ML1.m
function [Alpha, Beta, Gamma, Delta, s2u, s2v, suv] = AGGREGATE_ML1(Y, W,Z, Xagg)

% Purpose: 
% Analytic ML of aggregated covariate model (with instruments)
% -----------------------------------
% Model:
% Yti = Xti * Beta + Wti * Delta + uti
% Xti = Zti * Alpha  + Wti * Gamma + vti
% (uti, vti) ~ MVN[ (0,0), (su^2, suv; suv, sv^2) ]
% -----------------------------------
% Algorithm: 
% Qian (2010) 
% "Linear Regression Using Both Temporally Aggregated and Temporally Disaggregated Data: Revisited"
% -----------------------------------
% Usage:
% Y = Dependent variable (nT*1 vector
% W = Explantory variable W (nT*kw matrix
% Z = Explantory variable Z (nT*1 vector
% Xagg = Aggregated covariate X T*1 vector
% -----------------------------------
% Returns:
% Alpha = estimator of Alpha
% Beta = estimator of Beta
% Gamma = estimator of Gamma
% Delta = estimator of Delta
% s2u = estimator of su^2
% s2v = estimator of sv^2
% suv = estimator of suv
% 
% Written by Hang Qian, Iowa State University
% Contact me:  matlabist@gmail.com


% Analyze the data
[nobs,nregw] = size(W);
[nobs_test,nregz] = size(Z);
if nobs ~= nobs_test
    error('ErrorMismatched length of W and Z')
end

if nregz ~= 1
    error('Error: Analytic ML requires Zi be a scalar')
end

T = length(Xagg);
N = nobs / T;
if  N == 1
    error('Error: Aggregated Xagg should be shorter in length')
end
if N ~= floor(nobs / T)
    error('Error: Program only supports balanced aggregation')
end

% Construct aggregated variables
data = [Y,W,Z]; 
mat_temp = zeros(N,T);
data_agg = zeros(T,1+nregw+nregz);       
 for m = 1: (1+nregw+nregz)
    mat_temp(:) = data(:,m);
    data_agg(:,m) = sum(mat_temp);
 end
Yagg = data_agg(:,1);
Wagg = data_agg(:,2:nregw+1);
Zagg = data_agg(:,nregw+2:end);

% Reparameterize (A,B,C,D,E,F,G)OLS two equations
Reg1 = [Z,W];
OLS1 = (Reg1' * Reg1) \ (Reg1' * Y);
resid1 = Y - Reg1 * OLS1;
RSS1 = resid1' * resid1;
A = OLS1(1:nregz);
B = OLS1(nregz+1: nregz+nregw);
C = RSS1 / (N*T);

Reg2 = [Zagg,Wagg,Yagg];
OLS2 = (Reg2' * Reg2) \ (Reg2' * Xagg);
resid2 = Xagg - Reg2 * OLS2;
RSS2 = resid2' * resid2;
D = OLS2(1:nregz);
E = OLS2(nregz+1:nregz+nregw);
F = OLS2(nregz+nregw+1);
G = RSS2 / T;


% 7 equations solves for 7 unknowns
Alpha = D + A * F;
Beta = A / Alpha;
Gamma =  E + B * F;
Delta = (B * D - A * E) / Alpha;
s2u = (A^2 * G + N * C * D^2) / ( Alpha^2 * N);
s2v = C * F^2 + G / N;
suv = (N * C * D * F - A * G) / ( N * Alpha );


% Calculate standard error
OLS1_cov = inv(Reg1' * Reg1) * C;
OLS2_cov = inv(Reg2' * Reg2) * G;
V1 =blkdiag(OLS1_cov, 2/(N*T) * C^2);
V2 =blkdiag(OLS2_cov, 2/T * G^2);
V = blkdiag(V1,V2);

try
    disp('Try to calclulate analytic standard error (DELTA method)')    
    load Jacob7.mat;
    EE = E;
    Jacob_val = zeros(nregw*2+5);
    row = [1, 2, 3, 3+nregw, 3+nregw*2, 3+nregw*2+1, 3+nregw*2+2, 3+nregw*2+3];
    col = [1, 2, 2+nregw, 2+nregw+1,2+nregw+2, 4+nregw*2, 4+nregw*2+1, 4+nregw*2+2];
    for m = 1:7
        for n = 1:7            
             val = eval(Jacob(m,n)) ;
             if isscalar(val) && row(m+1)-row(m)>1 && col(n+1)-col(n)>1 
                 Jacob_val(row(m):row(m+1)-1, col(n):col(n+1)-1) = val * eye(nregw,nregw);
             else
                 Jacob_val(row(m):row(m+1)-1, col(n):col(n+1)-1) = val ;
             end
        end
    end
    
    cov_mat = Jacob_val * V *Jacob_val';
    theta_std=sqrt(diag(cov_mat));
    Alpha_std = theta_std(1:nregz);
    Beta_std = theta_std(nregz+1);
    Gamma_std = theta_std(nregz+2:nregz+nregw+1);
    Delta_std = theta_std(nregz+nregw+2:nregz+nregw*2+1);
    s2u_std = theta_std(nregz+nregw*2+2);
    s2v_std = theta_std(nregz+nregw*2+3);
    suv_std = theta_std(nregz+nregw*2+4);
    disp('Analytic standard error successfully calculated')
    
catch
    disp('Analytic standard error failed, switch to Monte Carlo standard error')
    ndraws = 3000;
    Alpha_sim = zeros(nregz,ndraws);
    Beta_sim = zeros(1,ndraws);
    Gamma_sim = zeros(nregw,ndraws);
    Delta_sim = zeros(nregw,ndraws);
    s2u_sim = zeros(1,ndraws);
    s2v_sim = zeros(1,ndraws);
    suv_sim = zeros(1,ndraws);
    
    P = chol(V)';
    para = [A; B; C; D; E; F; G];
    
    for r = 1:ndraws
        para_sim = para + P * randn(nregz+nregw*2+4,1);
        A = para_sim(1:nregz);
        B = para_sim(nregz+1: nregz+nregw);
        C = para_sim(nregz+nregw+1);
        D = para_sim(nregz+nregw+2);
        E = para_sim(nregz+nregw+3:nregz+nregw*2+2);
        F = para_sim(nregz+nregw*2+3);
        G = para_sim(nregz+nregw*2+4);
        
        common = D + A * F;
        Alpha_sim(:,r) = common;
        Beta_sim(:,r) = A / common;
        Gamma_sim(:,r) =  E + B * F;
        Delta_sim(:,r) = (B * D - A * E) / common;
        s2u_sim(:,r) = (A^2 * G + N * C * D^2) / ( common^2 * N);
        s2v_sim(:,r) = C * F^2 + G / N;
        suv_sim(:,r) = (N * C * D * F - A * G) / ( N * common );
    end
    
    Alpha_std = std(Alpha_sim,0,2);
    Beta_std = std(Beta_sim,0,2);
    Gamma_std = std(Gamma_sim,0,2);
    Delta_std = std(Delta_sim,0,2);
    s2u_std =std(s2u_sim,0,2);
    s2v_std = std(s2v_sim,0,2);
    suv_std = std(suv_sim,0,2);
end
    


% Display the estimation results
disp('Analytic MLE ')
disp('------Eq(1)------')
result=cell(nregw+3,3);
result(1,:)={'Variable',' Estimator',' SE'};  
result(2,1)={'Beta '};
result(2,2:3)={Beta,Beta_std};
for m=1:nregw
    result(m+2,1)={['Delta(', num2str(m), ')']};
    result(m+2,2:3)={Delta(m), Delta_std(m)};
end
result(end,1)={'su^2'};
result(end,2:3)={s2u,s2u_std};
disp(result)

disp(' ')
disp('------Eq(2)------')
result=cell(nregz+nregw+3,3);
result(1,:)={'Variable',' Estimator',' SE'};         
for m=1:nregz
    result(m+1,1)={['Alpha(',num2str(m),')']};
    result(m+1,2:3)={Alpha(m),Alpha_std(m)};
end
for m=1:nregw
    result(m+nregz+1,1)={['Gamma(',num2str(m),')']};
    result(m+nregz+1,2:3)={Gamma(m),Gamma_std(m)};
end
result(end-1,1)={'sv^2'};
result(end-1,2:3)={s2v,s2v_std};
result(end,1)={'suv'};
result(end,2:3)={suv,suv_std};
disp(result)

end

Contact us