% A Matlab/BACC file to illustrate estimation of a pair of seemingly
% unrelated regressions using (a) OLS, (b) GLS, (c) Bayesian estimation 
% with diffuse prior distributions.  
% The pair of equations to be estimated are factor demand equations
% derived from a generalized Leontief cost function.
%*******************************************************************
clear
% Generate artificial data. 
T=100;
setSeedConstant
q = 1+99*rand(T,1);  % uniform distribution on (1,100)
w1 = 1+99*rand(T,1);
w2 = 1+99*rand(T,1);
tempx = [q w1 w2];
corrcoef(tempx)
w2ow1 = w2./w1;  % elementwise division
w1ow2 = w1./w2;
tempx = [w1 w2 w2ow1 w1ow2];
corrcoef(tempx)
beta1 = 1;  % parameter values
beta2 = 1;
beta3 = 1;
y1det = beta1*q + beta2*q.*sqrt(w2ow1);
sigmasq1 = .5*var(y1det)
epsilon1 = sqrt(sigmasq1)*randn(T,1);
mean(epsilon1)
var(epsilon1)
tempe = sqrt(sigmasq1)*randn(T,1);
epsilon2 = .8*epsilon1 + .6*tempe;
mean(epsilon2)
var(epsilon2)
corrcoef(epsilon1,epsilon2)
qsrw2 = q.*sqrt(w2ow1);
qsrw1 = q.*sqrt(w1ow2);
tempx = [q(1:10,1) w2ow1(1:10,1) qsrw2(1:10,1) w1ow2(1:10,1) qsrw1(1:10,1)];
y1 = beta1*q + beta2*qsrw2 + epsilon1;
y2 = beta2*qsrw1 + beta3*q + epsilon2;
tempx = [q(1:10) qsrw2(1:10) epsilon1(1:10) y1(1:10)];
texpx = [qsrw1(1:10) q(1:10) epsilon2(1:10) y2(1:10)];
% OLS estimation
x1 = [q qsrw2];
x2 = [qsrw1 q];
a1 = (x1'*x1)\x1'*y1  % parameter estimates
a2 = (x2'*x2)\x2'*y2
yhat1 = x1*a1;
yhat2 = x2*a2;
resid1 = y1-yhat1;
resid2 = y2-yhat2;
s2y1 = resid1'*resid1/(T-2) % unbiased estimate of variance
s2y2 = resid2'*resid2/(T-2)
cvm1 = s2y1*inv(x1'*x1)
cvm2 = s2y2*inv(x2'*x2)
se1 = sqrt(diag(cvm1))  % standard errors
se2 = sqrt(diag(cvm2))
t1 = a1./se1  % t-statistics
t2 = a2./se2
% FGLS estimation
OLSres = [resid1 resid2];
SIGMAhat = (1/T)*OLSres'*OLSres % consistent estimate of cov. matrix
x = [x1 zeros(T,2); zeros(T,2) x2];
size(x)
OmegaInv = kron(inv(SIGMAhat),eye(T)); % 
OmegaInv(1:3, 1:3)
OmegaInv(1:3, T+1:T+3)
OmegaInv(2*T-2:2*T, 2*T-2:2*T)
covbetahat = inv(x'*OmegaInv*x)
se_fgls = sqrt(diag(covbetahat))
y = [y1;y2];
size(y)
betahat = covbetahat*x'*OmegaInv*y % fgls parameter estimates
tstat_fgls = betahat./se_fgls      % fgls t-statistics
% Bayesian estimation with cross-equation parameter restriction and
% diffuse prior distribution
z = [x1 zeros(T,1); zeros(T,1) x2]; % design matrix 
z(1:2,:)
z(2*T-1:2*T,:)
prioralpha = zeros(3,1);    % diffuse normal prior for slopes
priorstda =[100, 100, 100]';
priorvara = priorstda.^2;
priorCVMa = diag(priorvara);
priorHa = inv(priorCVMa);
nu_ = 4;    % diffuse Wishart prior for precision matrix of disturbances
s11_ = (1/3)*var(y1)
s22_ = (1/3)*var(y2)
S_ = [s11_ 0; 0 s22_] 
disp('Create model instance');
mi1 = minst('nlm','alpha','h',prioralpha,priorHa,nu_,S_,z,y); 
disp('Simulate model 1');
setSeedConstant
postsim(mi1,11000,1); % make 11000 draws from posterior distribution
filt = [1001:11000];  % discard burn-in draws
postfilter(mi1,filt); 
disp('Extract model instance to a Matlab structure');
sim1 = extract(mi1)
disp('Save model instance to a binary file');
miSave(mi1,'genLeontief_mi1.bin');
disp('Compute marginal likelihood')
[ml1,mlnse1] = mlike(mi1)
disp('Compute moments for parameters');
disp('alpha1')
expect1(sim1.logWeightPost,sim1.alpha(1,1,:))
disp('alpha2');
expect1(sim1.logWeightPost,sim1.alpha(2,1,:))
disp('alpha3');
expect1(sim1.logWeightPost,sim1.alpha(3,1,:))
disp('h1');
expect1(sim1.logWeightPost,sim1.h(1,1,:))
disp('h2');
expect1(sim1.logWeightPost,sim1.h(2,2,:))
disp('Plot sequence of draws from 1001 to 10000')
alpha1=squeeze(sim1.alpha(1,1,:))';
plot(alpha1)
xlabel('saved iteration','Fontsize',18);
ylabel('\alpha_1','Fontsize',20);
alpha2=squeeze(sim1.alpha(2,1,:))';
plot(alpha2)
xlabel('saved iteration','Fontsize',18);
ylabel('\alpha_2','Fontsize',20);
alpha3=squeeze(sim1.alpha(3,1,:))';
plot(alpha3);
xlabel('saved iteration','Fontsize',18);
ylabel('\alpha_3','Fontsize',20);
disp('Kernel smoothing and plotting of simulations')
[alpha1,pdf] = weightedSmooth(sim1.logWeightPost,sim1.alpha(1,1,:));
plot(alpha1,pdf)
xlabel('\alpha_1','Fontsize',20);
ylabel('P(\alpha_1|data)','Fontsize',20);
[alpha2,pdf] = weightedSmooth(sim1.logWeightPost,sim1.alpha(2,1,:));
plot(alpha2,pdf)
xlabel('\alpha_2','Fontsize',20);
ylabel('P(\alpha_2|data)','Fontsize',20);
[alpha3,pdf] = weightedSmooth(sim1.logWeightPost,sim1.alpha(3,1,:));
plot(alpha3,pdf)
xlabel('\alpha_3','Fontsize',20);
ylabel('P(\alpha_3|data)','Fontsize',20);
quit;