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_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]; 
    RSS_mat = resid_mat' * resid_mat;
    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_mean_adjust = X_post_mean_mat + gap(ones(N,1),:);
    X_post_var_adjust = X_post_var * (eye(N-1) - 1/N);  

    Xsim_mat(1:end-1,:) = X_post_mean_adjust(1:end-1,:) + chol(X_post_var_adjust)' * randn(N-1,T);
    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




Contact us