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_GIBBS1.m
```function [Alpha, Beta, Gamma, Delta, s2u, s2v, suv] = AGGREGATE_GIBBS1(Y, W, Z, Xagg, ndraws, burn_in, initial_val)

% Purpose:
% Bayesian estimation of aggregated covariate model via Gibbs sampler
% -----------------------------------
% 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
% ndraws = number of draws in Gibbs sampler
% burn_in = burn in periods in Gibbs sampler
% -----------------------------------
% 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

if nargin < 5;  ndraws = 10000;   end
if nargin < 6;  burn_in = ndraws * 0.5;   end
if nargin < 7; initial_flag = 0; else initial_flag = 1; end

% Analyze the data
[nobs,nregw] = size(W);
[nobs_test,nregz] = size(Z);
nreg_eq1 = nregw+1;
nreg_eq2 = nregz + nregw;
npara = nreg_eq1 + nreg_eq2;

if nobs ~= nobs_test
error('ErrorMismatched length of W and Z')
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

%----------------------------------------
% Priors specification
%  Beta ~ N(mu,V)
%  (su^2, suv; suv, sv^2) ~ Wishart(Omega,df)
prior_mu = zeros(npara,1);
prior_V = 100*eye(npara);
prior_omega = eye(2);
prior_df = 1;

%-----------------------------------------

% Allocating memory and set initial values
Beta = zeros(1 , ndraws - burn_in);
Delta = zeros(nregw , ndraws - burn_in);
Alpha = zeros(nregz , ndraws - burn_in);
Gamma = zeros(nregw , ndraws - burn_in);
Sigma = zeros(2,2,ndraws-burn_in);
X_post_mean_mat = zeros(N,T);
X_SUR_sd=zeros(nobs*2, npara);

if initial_flag == 1
if length(initial_val) ~= nregz+nregw*2+4
error('Error: Loaded inital value has wrong length')
end
Alpha_use = initial_val(1:nregz);
Beta_use = initial_val(nregz+1);
Gamma_use = initial_val(nregz+2:nregz+nregw+1);
Delta_use = initial_val(nregz+nregw+2:nregz+nregw*2+1);
s2u_use = initial_val(nregz+nregw*2+2);
s2v_use = initial_val(nregz+nregw*2+3);
suv_use = initial_val(nregz+nregw*2+4);
Sigma_use = [ s2u_use, suv_use;suv_use,s2v_use];
inv_Sigma = inv(Sigma_use);
else
Betabig = prior_mu + chol(prior_V)' * randn(npara,1);
Beta_use = Betabig(1);
Delta_use = Betabig(2:nreg_eq1);
Alpha_use = Betabig(nreg_eq1+1 : nreg_eq1+nregz);
Gamma_use = Betabig(nreg_eq1+nregz+1:end);
inv_Sigma = WISHART_RND1(prior_omega, prior_df);
end

Xavg = Xagg' / N;
Xsim_mat = Xavg(ones(N,1),:);
Xsim = Xsim_mat(:);

XW = [Xsim,W];
ZW = [Z,W];

% Gibbs sampler starts
%H_status_bar = waitbar(0,'Simulation in progress   0%');
%percentage = floor(ndraws/100);
for r = 1:ndraws

%     if mod(r,percentage)==0
%        waitbar(r/nrep,H_status_bar,['Simulation in progress   ',num2str(floor(r/nrep*100)),'%']);
%     end

% Standard SUR model,sample from posterior
fit_eq1 =  Xsim * Beta_use + W * Delta_use;
fit_eq2 = Z *  Alpha_use + W * Gamma_use;
resid_mat = [ Y - fit_eq1 , Xsim - fit_eq2];
post_omega = inv( inv(prior_omega) + RSS_mat);
inv_Sigma = WISHART_RND1(post_omega, nobs + prior_df);
Sigma_use = inv(inv_Sigma);

% % Standard SUR model,sample from posterior (;;,)
% You can also use the following codes, but needs patience
% D = inv(prior_V);    d = prior_V \ prior_mu;
% for m = 1: nobs
%    Xbig = blkdiag( [Xsim(m), W(m,:)],  [Z(m,:), W(m,:)] );
%    D = D + Xbig' * inv_Sigma * Xbig;
%    d = d + Xbig' * inv_Sigma * [Y(m);Xsim(m)] ;
% end
% D = inv(D);
% post_mean = D * d;
% post_var = D;

P = chol(inv_Sigma);
Y_mat_sd = P *  [Y,Xsim]';
Y_SUR_sd = Y_mat_sd(:);
XW(:,1) = Xsim;
X_SUR_sd(1:2:end-1,1:nreg_eq1) = P(1,1) * XW;
X_SUR_sd(1:2:end-1,nreg_eq1+1:end) = P(1,2) * ZW;
X_SUR_sd(2:2:end  ,nreg_eq1+1:end) = P(2,2) * ZW;

D = inv( X_SUR_sd' * X_SUR_sd +  inv(prior_V) );
d = X_SUR_sd' * Y_SUR_sd + prior_V \ prior_mu;
post_mean = D * d;
post_var = D;
Betabig = post_mean + chol(post_var)' * randn(npara,1);
Beta_use = Betabig(1);
Delta_use = Betabig(2:nreg_eq1);
Alpha_use = Betabig(nreg_eq1+1 : nreg_eq1+nregz);
Gamma_use = Betabig(nreg_eq1+nregz+1:end);

%----------------
% Sample from the latent X
%----------------
cov11 = Sigma_use(2,2);
cov12 = Beta_use * Sigma_use(2,2) + Sigma_use(1,2);
cov22 = Beta_use^2 * Sigma_use(2,2) + Sigma_use(1,1) + 2 * Beta_use * Sigma_use(1,2);
fit_eq2 = Z *  Alpha_use + W * Gamma_use;
X_post_mean = fit_eq2 + cov12/cov22 * (Y - fit_eq2 * Beta_use - W * Delta_use);
X_post_var = cov11 - cov12^2 / cov22;

X_post_mean_mat(:) =  X_post_mean;
gap = (Xagg' - sum(X_post_mean_mat) ) / N;

X_post_var_adjust = X_post_var * (eye(N-1) - 1/N);

Xsim_mat(end,:) = Xagg' - sum( Xsim_mat(1:end-1,:));
Xsim = Xsim_mat(:);

% Record the draws
if r > burn_in
Beta(:, r - burn_in) = Beta_use;
Delta(:, r - burn_in) = Delta_use;
Alpha(:, r - burn_in) = Alpha_use;
Gamma(:, r - burn_in) = Gamma_use;
Sigma(:, :, r - burn_in) = Sigma_use;
end

end

% Compute mean and standard deviation
s2u = zeros(1, ndraws - burn_in);
s2v = zeros(1, ndraws - burn_in);
suv = zeros(1, ndraws - burn_in);
s2u(:) = Sigma(1,1,:);
s2v(:) = Sigma(2,2,:);
suv(:) = Sigma(1,2,:);
Beta_mean = mean(Beta,2);
Beta_std = std(Beta,0,2);
Delta_mean = mean(Delta,2);
Delta_std = std(Delta,0,2);
Alpha_mean = mean(Alpha,2);
Alpha_std = std(Alpha,0,2);
Gamma_mean = mean(Gamma,2);
Gamma_std = std(Gamma,0,2);

s2u_mean = mean(s2u,2);
s2u_std = std(s2u,0,2);
s2v_mean = mean(s2v,2);
s2v_std = std(s2v,0,2);
suv_mean = mean(suv,2);
suv_std = std(suv,0,2);

% Display the results
disp(' ')
disp('MCMC esitmation results ')
disp('------Eq(1)------')
result=cell(nregw+3,3);
result(1,:)={'Variable',' Estimator',' SE'};
result(2,1)={'Beta '};
result(2,2:3)={Beta_mean,Beta_std};
for m=1:nregw
result(m+2,1)={['Delta(',num2str(m),')']};
result(m+2,2:3)={Delta_mean(m),Delta_std(m)};
end
result(nregw+3,1)={'su^2'};
result(nregw+3,2:3)={s2u_mean,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_mean(m),Alpha_std(m)};
end
for m=1:nregw
result(nregz+m+1,1)={['Gamma(',num2str(m),')']};
result(nregz+m+1,2:3)={Gamma_mean(m),Gamma_std(m)};
end
result(nregz+nregw+2,1)={'sv^2'};
result(nregz+nregw+2,2:3)={s2v_mean,s2v_std};
result(nregz+nregw+3,1)={'suv'};
result(nregz+nregw+3,2:3)={suv_mean,suv_std};
disp(result)

%close(H_status_bar);

end

```