Code covered by the BSD License

Toolkit on Econometrics and Economics Teaching

Hang Qian (view profile)

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;
A = OLS1(1:nregz);
B = OLS1(nregz+1: nregz+nregw);

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

% 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)')
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

```