%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 1/ Xij: Trade flows from exporter i to importer j
% 2/ Xij13: CAN-US Trade Flow*1.3
% 3/ Xij23: CAN-US Trade Flow*1.6
% 4/ Yi: GDP of exporter
% 5/ Yj: GDP of importeur
% 6/ dijAvW: bilateral distance with Anderson and van Wincoop (2003)
% intra-distance
% 7/ dij13: bilateral distance with surface type intra-distance (using 1/3)
% 8/ dij23: bilateral distance with surface type intra-distance (using 2/3)
% 9/ popi: population of exporter
% 10/ wage: hourly worked of exporter
% 11/ vm_cancan = 1 if i=CAN and j=CAN (McCallum, 1995)
% 12/ vm_usus = 1 if i=US and j=US (McCallum, 1995)
% 13/ vm_canus = 1 if i=CAN and j=US (Anderson and van Wincoop, 2003)
% 14/ vm_uscan = 1 if i=US and j=CAN (Anderson and van Wincoop, 2003)
% 15/ vm0 = 1 if trade flow=0
% 16/ dcan = 1 if i=CAN
% 17/ dus = 1 if i=US
% 18-95/ Feenstra (2005) fixed effects
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Data
%
clear all;
load data.txt;
lnXij = log(data(:,1))-log(1000000);
lnXij13 = log(data(:,2))-log(1000000);
lnXij16 = log(data(:,3))-log(1000000);
lnYi = log(data(:,4));
lnYj = log(data(:,5));
lndijAvW = log(data(:,6));
lndij13 = log(data(:,7));
lndij23 = log(data(:,8));
lnpopi = log(data(:,9));
lnwage = log(data(:,10));
CAN = data(:,11);
US = data(:,12);
vmCAN = data(:,13);
vmUS = data(:,14);
vm0 = data(:,15);
n2 = length(data);
n = sqrt(n2);
info.lflag = 0;
%
% Interaction matrix
%
P = zeros(n,1);
for i=1:n;
P(i,1) = data(i*n,9);
end;
Sp = sum(P);
Ln = (1/Sp)*P;
Wt = kron((ones(n,n))*(diag(Ln)),eye(n,n));
Wd = diag(diag(Wt));
W = Wt - Wd;
%
lnLi = lnpopi - log(Sp)*ones(n2,1);
%
% Endogeneous variable
%
y = lnXij - lnYi - lnYj;
%
% McCallum
%
x1 = [ones(n2,1) lnYi lnYj lndijAvW CAN US vm0];
vnames = strvcat('lnXij','constant','lnYi','lnYj','lndij','VMCAN','VMUS','vm0');
res1 = hwhite(lnXij,x1);
prt(res1,vnames);
res1.beta./res1.tstat
bCAN1 = exp(res1.beta(5))
bUS1 = exp(res1.beta(6))
bav1 = sqrt(bCAN1*bUS1)
resmor1 = moran(lnXij,x1,W);
prt(resmor1);
J1 = [exp(res1.beta(5)) 0
       0 exp(res1.beta(6))
      ((exp(res1.beta(5)+res1.beta(6)))/(2*sqrt(exp(res1.beta(5)+res1.beta(6))))) ((exp(res1.beta(5)+res1.beta(6)))/(2*sqrt(exp(res1.beta(5)+res1.beta(6)))))];
V1 = res1.minfo;
var1 = J1*V1(5:6,5:6)*J1';
sigbCAN1 = sqrt(var1(1,1))
sigbUS1 = sqrt(var1(2,2))
sigbav1 = sqrt(var1(3,3))
tbCAN1 = abs(bCAN1 / sigbCAN1);
tbUS1 = abs(bUS1 / sigbUS1);
tbav1 = abs(bav1 / sigbav1);
pbCAN1 = 2*(1 - norm_cdf (tbCAN1))
pbUS1 = 2*(1 - norm_cdf (tbUS1))
pbav1 = 2*(1 - norm_cdf (tbav1))
%
x2 = [ones(n2,1) lnYi lnYj lndij13 CAN US vm0];
vnames = strvcat('lnXij','constant','lnYi','lnYj','lndij','VMCAN','VMUS','vm0');
res2 = hwhite(lnXij,x2);
prt(res2,vnames);
res2.beta./res2.tstat
bCAN2 = exp(res2.beta(5))
bUS2 = exp(res2.beta(6))
bav2 = sqrt(bCAN2*bUS2)
resmor2 = moran(lnXij,x2,W);
prt(resmor2);
J2 = [exp(res2.beta(5)) 0
       0 exp(res2.beta(6))
      ((exp(res2.beta(5)+res2.beta(6)))/(2*sqrt(exp(res2.beta(5)+res2.beta(6))))) ((exp(res2.beta(5)+res2.beta(6)))/(2*sqrt(exp(res2.beta(5)+res2.beta(6)))))];
V2 = res2.minfo;
var2 = J2*V2(5:6,5:6)*J2';
sigbCAN2 = sqrt(var2(1,1))
sigbUS2 = sqrt(var2(2,2))
sigbav2 = sqrt(var2(3,3))
tbCAN2 = abs(bCAN2 / sigbCAN2);
tbUS2 = abs(bUS2 / sigbUS2);
tbav2 = abs(bav2 / sigbav2);
pbCAN2 = 2*(1 - norm_cdf (tbCAN2))
pbUS2 = 2*(1 - norm_cdf (tbUS2))
pbav2 = 2*(1 - norm_cdf (tbav2))
%
x3 = [ones(n2,1) lnYi lnYj lndij23 CAN US vm0];
vnames = strvcat('lnXij','constant','lnYi','lnYj','lndij','VMCAN','VMUS','vm0');
res3 = hwhite(lnXij,x3);
prt(res3,vnames);
res3.beta./res3.tstat
bCAN3 = exp(res3.beta(5))
bUS3 = exp(res3.beta(6))
bav3 = sqrt(bCAN3*bUS3)
resmor3 = moran(lnXij,x3,W);
prt(resmor3);
J3 = [exp(res3.beta(5)) 0
       0 exp(res3.beta(6))
      (exp(res3.beta(5)+res3.beta(6)))/(2*sqrt(exp(res3.beta(5)+res3.beta(6)))) ((exp(res3.beta(5)+res3.beta(6)))/(2*sqrt(exp(res3.beta(5)+res3.beta(6)))))];
V3 = res3.minfo;
var3 = J3*V3(5:6,5:6)*J3';
sigbCAN3 = sqrt(var3(1,1))
sigbUS3 = sqrt(var3(2,2))
sigbav3 = sqrt(var3(3,3))
tbCAN3 = abs(bCAN3 / sigbCAN3);
tbUS3 = abs(bUS3 / sigbUS3);
tbav3 = abs(bav3 / sigbav3);
pbCAN3 = 2*(1 - norm_cdf (tbCAN3))
pbUS3 = 2*(1 - norm_cdf (tbUS3))
pbav3 = 2*(1 - norm_cdf (tbav3))
%
% OLS Feenstra
%
xf1 = [ones(n2,1) lndijAvW vmCAN+vmUS vm0 data(:,18:95)];
resf1 = hwhite(y,xf1);
prt(resf1);
resf1.beta./resf1.tstat
resfmor1 = moran(y,xf1,W);
prt(resfmor1);
bavf1 = exp(-resf1.beta(3))
Jf1 = -(exp(-resf1.beta(3)))/(2*sqrt(exp(-resf1.beta(3))));
Vf1 = resf1.minfo;
varf1 = Jf1*Vf1(3,3)*Jf1;
sigbavf1 = sqrt(varf1)
tbavf1 = abs(bavf1 / sigbavf1);
pbavf1 = 2*(1 - norm_cdf (tbavf1))
%
xf2 = [ones(n2,1) lndij13 vmCAN+vmUS vm0 data(:,18:95)];
resf2 = hwhite(y,xf2);
prt(resf2);
resf2.beta./resf2.tstat
resfmor2 = moran(y,xf2,W);
prt(resfmor2);
bavf2 = exp(-resf2.beta(3))
Jf2 = -(exp(-resf2.beta(3)))/(2*sqrt(exp(-resf2.beta(3))));
Vf2 = resf2.minfo;
varf2 = Jf2*Vf2(3,3)*Jf2;
sigbavf2 = sqrt(varf2)
tbavf2 = abs(bavf2 / sigbavf2);
pbavf2 = 2*(1 - norm_cdf (tbavf2))
%
xf3 = [ones(n2,1) lndij23 vmCAN+vmUS vm0 data(:,18:95)];
resf3 = hwhite(y,xf3);
prt(resf3);
resf3.beta./resf3.tstat
resfmor3 = moran(y,xf3,W);
prt(resfmor3);
bavf3 = exp(-resf3.beta(3))
Jf3 = -(exp(-resf3.beta(3)))/(2*sqrt(exp(-resf3.beta(3))));
Vf3 = resf3.minfo;
varf3 = Jf3*Vf3(3,3)*Jf3;
sigbavf3 = sqrt(varf3)
tbavf3 = abs(bavf3 / sigbavf3);
pbavf3 = 2*(1 - norm_cdf (tbavf3))
%
% SAR model
%
xsar1 = [ones(n2,1) (eye(n2,n2) - Wt)*lndijAvW lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressar1 = sargrav(y,xsar1,W,info);
prt(ressar1,vnames);
[ressar1.beta
 ressar1.rho]./ressar1.tstat
% Criteria information
logliksar1 = ressar1.lik;
AICsar1 = logliksar1*(2/n2) + ressar1.nvar*(2/n2)
SCsar1 = logliksar1*(2/n2) + ressar1.nvar*(log(n2)/n2)
% Borders
bsar1 = exp(inv(eye(n2,n2) - ressar1.rho*Wt)*(ressar1.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarCANCAN1 = bsar1(1)
bsarCANUS1 = bsar1(11)
bsarUSUS1 = bsar1(411)
bsarUSCAN1 = bsar1(401)
bsarCAN1 = bsarCANCAN1/bsarCANUS1
bsarUS1 = bsarUSUS1/bsarUSCAN1
bsarav1 = sqrt(bsarCAN1*bsarUS1)
% p-value of partial border
Jsar1 = [(inv(eye(n2,n2) - ressar1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressar1.beta(4)*(inv(eye(n2,n2) - ressar1.rho*Wt))*Wt*(inv(eye(n2,n2) - ressar1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsar1 = [ressar1.minfo(4,4) ressar1.minfo(4,6)
    ressar1.minfo(6,4) ressar1.minfo(6,6)];
varsar1 = Jsar1*Vsar1*Jsar1';
sigbsarCANCAN1 = sqrt(varsar1(1,1))
sigbsarCANUS1 = sqrt(varsar1(11,11))
sigbsarUSUS1 = sqrt(varsar1(411,411))
sigbsarUSCAN1 = sqrt(varsar1(401,401))
tbsarCANCAN1 = abs(bsarCANCAN1 / sigbsarCANCAN1);
tbsarCANUS1 = abs(bsarCANUS1 / sigbsarCANUS1);
tbsarUSUS1 = abs(bsarUSUS1 / sigbsarUSUS1);
tbsarUSCAN1 = abs(bsarUSCAN1 / sigbsarUSCAN1);
pbsarCANCAN1 = 2*(1 - norm_cdf (tbsarCANCAN1))
pbsarCANUS1 = 2*(1 - norm_cdf (tbsarCANUS1))
pbsarUSUS1 = 2*(1 - norm_cdf (tbsarUSUS1))
pbsarUSCAN1 = 2*(1 - norm_cdf (tbsarUSCAN1))
% p-value of full border
Jsar1f = [1/bsarCANUS1  -(bsarCANCAN1)/((bsarCANUS1)^2)  0  0
    0  0  1/bsarUSCAN1  -(bsarUSUS1)/((bsarUSCAN1)^2)
    (1/(2*sqrt((bsarCANCAN1*bsarUSUS1)/(bsarCANUS1*bsarUSCAN1))))*((bsarUSUS1)/(bsarCANUS1*bsarUSCAN1))  -(1/(2*sqrt((bsarCANCAN1*bsarUSUS1)/(bsarCANUS1*bsarUSCAN1))))*((bsarCANCAN1*bsarUSUS1*bsarUSCAN1)/((bsarCANUS1*bsarUSCAN1)^2))  -(1/(2*sqrt((bsarCANCAN1*bsarUSUS1)/(bsarCANUS1*bsarUSCAN1))))*((bsarCANCAN1*bsarUSUS1*bsarCANUS1)/((bsarCANUS1*bsarUSCAN1)^2))  (1/(2*sqrt((bsarCANCAN1*bsarUSUS1)/(bsarCANUS1*bsarUSCAN1))))*((bsarCANCAN1)/(bsarCANUS1*bsarUSCAN1))];
subvarsar1 = [varsar1(1,1) varsar1(1,11) varsar1(1,411) varsar1(1,401)
              varsar1(11,1) varsar1(11,11) varsar1(11,411) varsar1(11,401)
              varsar1(411,1) varsar1(411,11) varsar1(411,411) varsar1(411,401)
              varsar1(401,1) varsar1(401,11) varsar1(401,411) varsar1(401,401)];
varsar1f = Jsar1f*subvarsar1*Jsar1f';
sigbsarCAN1 = sqrt(varsar1f(1,1))
sigbsarUS1 = sqrt(varsar1f(2,2))
sigbsarav1 = sqrt(varsar1f(3,3))
tbsarCAN1 = abs(bsarCAN1 / sigbsarCAN1);
tbsarUS1 = abs(bsarUS1 / sigbsarUS1);
tbsarav1 = abs(bsarav1 / sigbsarav1);
pbsarCAN1 = 2*(1 - norm_cdf (tbsarCAN1))
pbsarUS1 = 2*(1 - norm_cdf (tbsarUS1))
pbsarav1 = 2*(1 - norm_cdf (tbsarav1))
%
%
xsar2 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij13 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressar2 = sargrav(y,xsar2,W,info);
prt(ressar2,vnames);
[ressar2.beta
 ressar2.rho]./ressar2.tstat
% Criteria information
logliksar2 = ressar2.lik;
AICsar2 = logliksar2*(2/n2) + ressar2.nvar*(2/n2)
SCsar2 = logliksar2*(2/n2) + ressar2.nvar*(log(n2)/n2)
% Borders
bsar2 = exp(inv(eye(n2,n2) - ressar2.rho*Wt)*(ressar2.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarCANCAN2 = bsar2(1)
bsarCANUS2 = bsar2(11)
bsarUSUS2 = bsar2(411)
bsarUSCAN2 = bsar2(401)
bsarCAN2 = bsarCANCAN2/bsarCANUS2
bsarUS2 = bsarUSUS2/bsarUSCAN2
bsarav2 = sqrt(bsarCAN2*bsarUS2)
% p-value of partial border
Jsar2 = [(inv(eye(n2,n2) - ressar2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressar2.beta(4)*(inv(eye(n2,n2) - ressar2.rho*Wt))*Wt*(inv(eye(n2,n2) - ressar2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsar2 = [ressar2.minfo(4,4) ressar2.minfo(4,6)
    ressar2.minfo(6,4) ressar2.minfo(6,6)];
varsar2 = Jsar2*Vsar2*Jsar2';
sigbsarCANCAN2 = sqrt(varsar2(1,1))
sigbsarCANUS2 = sqrt(varsar2(11,11))
sigbsarUSUS2 = sqrt(varsar2(411,411))
sigbsarUSCAN2 = sqrt(varsar2(401,401))
tbsarCANCAN2 = abs(bsarCANCAN2 / sigbsarCANCAN2);
tbsarCANUS2 = abs(bsarCANUS2 / sigbsarCANUS2);
tbsarUSUS2 = abs(bsarUSUS2 / sigbsarUSUS2);
tbsarUSCAN2 = abs(bsarUSCAN2 / sigbsarUSCAN2);
pbsarCANCAN2 = 2*(1 - norm_cdf (tbsarCANCAN2))
pbsarCANUS2 = 2*(1 - norm_cdf (tbsarCANUS2))
pbsarUSUS2 = 2*(1 - norm_cdf (tbsarUSUS2))
pbsarUSCAN2 = 2*(1 - norm_cdf (tbsarUSCAN2))
% p-value of full border
Jsar2f = [1/bsarCANUS2  -(bsarCANCAN2)/((bsarCANUS2)^2)  0  0
    0  0  1/bsarUSCAN2  -(bsarUSUS2)/((bsarUSCAN2)^2)
    (1/(2*sqrt((bsarCANCAN2*bsarUSUS2)/(bsarCANUS2*bsarUSCAN2))))*((bsarUSUS2)/(bsarCANUS2*bsarUSCAN2))  -(1/(2*sqrt((bsarCANCAN2*bsarUSUS2)/(bsarCANUS2*bsarUSCAN2))))*((bsarCANCAN2*bsarUSUS2*bsarUSCAN2)/((bsarCANUS2*bsarUSCAN2)^2))  -(1/(2*sqrt((bsarCANCAN2*bsarUSUS2)/(bsarCANUS2*bsarUSCAN2))))*((bsarCANCAN2*bsarUSUS2*bsarCANUS2)/((bsarCANUS2*bsarUSCAN2)^2))  (1/(2*sqrt((bsarCANCAN2*bsarUSUS2)/(bsarCANUS2*bsarUSCAN2))))*((bsarCANCAN2)/(bsarCANUS2*bsarUSCAN2))];
subvarsar2 = [varsar2(1,1) varsar2(1,11) varsar2(1,411) varsar2(1,401)
              varsar2(11,1) varsar2(11,11) varsar2(11,411) varsar2(11,401)
              varsar2(411,1) varsar2(411,11) varsar2(411,411) varsar2(411,401)
              varsar2(401,1) varsar2(401,11) varsar2(401,411) varsar2(401,401)];
varsar2f = Jsar2f*subvarsar2*Jsar2f';
sigbsarCAN2 = sqrt(varsar2f(1,1))
sigbsarUS2 = sqrt(varsar2f(2,2))
sigbsarav2 = sqrt(varsar2f(3,3))
tbsarCAN2 = abs(bsarCAN2 / sigbsarCAN2);
tbsarUS2 = abs(bsarUS2 / sigbsarUS2);
tbsarav2 = abs(bsarav2 / sigbsarav2);
pbsarCAN2 = 2*(1 - norm_cdf (tbsarCAN2))
pbsarUS2 = 2*(1 - norm_cdf (tbsarUS2))
pbsarav2 = 2*(1 - norm_cdf (tbsarav2))
%
%
xsar3 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij23 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressar3 = sargrav(y,xsar3,W,info);
prt(ressar3,vnames);
[ressar3.beta
 ressar3.rho]./ressar3.tstat
% Criteria information
logliksar3 = ressar3.lik;
AICsar3 = logliksar3*(2/n2) + ressar3.nvar*(2/n2)
SCsar3 = logliksar3*(2/n2) + ressar3.nvar*(log(n2)/n2)
% Borders
bsar3 = exp(inv(eye(n2,n2) - ressar3.rho*Wt)*(ressar3.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarCANCAN3 = bsar3(1)
bsarCANUS3 = bsar3(11)
bsarUSUS3 = bsar3(411)
bsarUSCAN3 = bsar3(401)
bsarCAN3 = bsarCANCAN3/bsarCANUS3
bsarUS3 = bsarUSUS3/bsarUSCAN3
bsarav3 = sqrt(bsarCAN3*bsarUS3)
% p-value of partial border
Jsar3 = [(inv(eye(n2,n2) - ressar3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressar3.beta(4)*(inv(eye(n2,n2) - ressar3.rho*Wt))*Wt*(inv(eye(n2,n2) - ressar3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsar3 = [ressar3.minfo(4,4) ressar3.minfo(4,6)
    ressar3.minfo(6,4) ressar3.minfo(6,6)];
varsar3 = Jsar3*Vsar3*Jsar3';
sigbsarCANCAN3 = sqrt(varsar3(1,1))
sigbsarCANUS3 = sqrt(varsar3(11,11))
sigbsarUSUS3 = sqrt(varsar3(411,411))
sigbsarUSCAN3 = sqrt(varsar3(401,401))
tbsarCANCAN3 = abs(bsarCANCAN3 / sigbsarCANCAN3);
tbsarCANUS3 = abs(bsarCANUS3 / sigbsarCANUS3);
tbsarUSUS3 = abs(bsarUSUS3 / sigbsarUSUS3);
tbsarUSCAN3 = abs(bsarUSCAN3 / sigbsarUSCAN3);
pbsarCANCAN3 = 2*(1 - norm_cdf (tbsarCANCAN3))
pbsarCANUS3 = 2*(1 - norm_cdf (tbsarCANUS3))
pbsarUSUS3 = 2*(1 - norm_cdf (tbsarUSUS3))
pbsarUSCAN3 = 2*(1 - norm_cdf (tbsarUSCAN3))
% p-value of full border
Jsar3f = [1/bsarCANUS3  -(bsarCANCAN3)/((bsarCANUS3)^2)  0  0
    0  0  1/bsarUSCAN3  -(bsarUSUS3)/((bsarUSCAN3)^2)
    (1/(2*sqrt((bsarCANCAN3*bsarUSUS3)/(bsarCANUS3*bsarUSCAN3))))*((bsarUSUS3)/(bsarCANUS3*bsarUSCAN3))  -(1/(2*sqrt((bsarCANCAN3*bsarUSUS3)/(bsarCANUS3*bsarUSCAN3))))*((bsarCANCAN3*bsarUSUS3*bsarUSCAN3)/((bsarCANUS3*bsarUSCAN3)^2))  -(1/(2*sqrt((bsarCANCAN3*bsarUSUS3)/(bsarCANUS3*bsarUSCAN3))))*((bsarCANCAN3*bsarUSUS3*bsarCANUS3)/((bsarCANUS3*bsarUSCAN3)^2))  (1/(2*sqrt((bsarCANCAN3*bsarUSUS3)/(bsarCANUS3*bsarUSCAN3))))*((bsarCANCAN3)/(bsarCANUS3*bsarUSCAN3))];
subvarsar3 = [varsar3(1,1) varsar3(1,11) varsar3(1,411) varsar3(1,401)
              varsar3(11,1) varsar3(11,11) varsar3(11,411) varsar3(11,401)
              varsar3(411,1) varsar3(411,11) varsar3(411,411) varsar3(411,401)
              varsar3(401,1) varsar3(401,11) varsar3(401,411) varsar3(401,401)];
varsar3f = Jsar3f*subvarsar3*Jsar3f';
sigbsarCAN3 = sqrt(varsar3f(1,1))
sigbsarUS3 = sqrt(varsar3f(2,2))
sigbsarav3 = sqrt(varsar3f(3,3))
tbsarCAN3 = abs(bsarCAN3 / sigbsarCAN3);
tbsarUS3 = abs(bsarUS3 / sigbsarUS3);
tbsarav3 = abs(bsarav3 / sigbsarav3);
pbsarCAN3 = 2*(1 - norm_cdf (tbsarCAN3))
pbsarUS3 = 2*(1 - norm_cdf (tbsarUS3))
pbsarav3 = 2*(1 - norm_cdf (tbsarav3))
%
% GSM model
%
xsac1 = [ones(n2,1) (eye(n2,n2) - Wt)*lndijAvW lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressac1 = sacgrav(y,xsac1,W,W,info);
prt(ressac1,vnames);
[ressac1.beta
 ressac1.rho
 ressac1.lam]./ressac1.tstat
% Criteria information
logliksac1 = ressac1.lik;
AICsac1 = logliksac1*(2/n2) + ressac1.nvar*(2/n2)
SCsac1 = logliksac1*(2/n2) + ressac1.nvar*(log(n2)/n2)
% Borders
bsac1 = exp(inv(eye(n2,n2) - ressac1.rho*Wt)*(ressac1.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsacCANCAN1 = bsac1(1)
bsacCANUS1 = bsac1(11)
bsacUSUS1 = bsac1(411)
bsacUSCAN1 = bsac1(401)
bsacCAN1 = bsacCANCAN1/bsacCANUS1
bsacUS1 = bsacUSUS1/bsacUSCAN1
bsacav1 = sqrt(bsacCAN1*bsacUS1)
% p-value of partial border
Jsac1 = [(inv(eye(n2,n2) - ressac1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressac1.beta(4)*(inv(eye(n2,n2) - ressac1.rho*Wt))*Wt*(inv(eye(n2,n2) - ressac1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsac1 = [ressac1.minfo(4,4) ressac1.minfo(4,6)
    ressac1.minfo(6,4) ressac1.minfo(6,6)];
varsac1 = Jsac1*Vsac1*Jsac1';
sigbsacCANCAN1 = sqrt(varsac1(1,1))
sigbsacCANUS1 = sqrt(varsac1(11,11))
sigbsacUSUS1 = sqrt(varsac1(411,411))
sigbsacUSCAN1 = sqrt(varsac1(401,401))
tbsacCANCAN1 = abs(bsacCANCAN1 / sigbsacCANCAN1);
tbsacCANUS1 = abs(bsacCANUS1 / sigbsacCANUS1);
tbsacUSUS1 = abs(bsacUSUS1 / sigbsacUSUS1);
tbsacUSCAN1 = abs(bsacUSCAN1 / sigbsacUSCAN1);
pbsacCANCAN1 = 2*(1 - norm_cdf (tbsacCANCAN1))
pbsacCANUS1 = 2*(1 - norm_cdf (tbsacCANUS1))
pbsacUSUS1 = 2*(1 - norm_cdf (tbsacUSUS1))
pbsacUSCAN1 = 2*(1 - norm_cdf (tbsacUSCAN1))
% p-value of full border
Jsac1f = [1/bsacCANUS1  -(bsacCANCAN1)/((bsacCANUS1)^2)  0  0
    0  0  1/bsacUSCAN1  -(bsacUSUS1)/((bsacUSCAN1)^2)
    (1/(2*sqrt((bsacCANCAN1*bsacUSUS1)/(bsacCANUS1*bsacUSCAN1))))*((bsacUSUS1)/(bsacCANUS1*bsacUSCAN1))  -(1/(2*sqrt((bsacCANCAN1*bsacUSUS1)/(bsacCANUS1*bsacUSCAN1))))*((bsacCANCAN1*bsacUSUS1*bsacUSCAN1)/((bsacCANUS1*bsacUSCAN1)^2))  -(1/(2*sqrt((bsacCANCAN1*bsacUSUS1)/(bsacCANUS1*bsacUSCAN1))))*((bsacCANCAN1*bsacUSUS1*bsacCANUS1)/((bsacCANUS1*bsacUSCAN1)^2))  (1/(2*sqrt((bsacCANCAN1*bsacUSUS1)/(bsacCANUS1*bsacUSCAN1))))*((bsacCANCAN1)/(bsacCANUS1*bsacUSCAN1))];
subvarsac1 = [varsac1(1,1) varsac1(1,11) varsac1(1,411) varsac1(1,401)
              varsac1(11,1) varsac1(11,11) varsac1(11,411) varsac1(11,401)
              varsac1(411,1) varsac1(411,11) varsac1(411,411) varsac1(411,401)
              varsac1(401,1) varsac1(401,11) varsac1(401,411) varsac1(401,401)];
varsac1f = Jsac1f*subvarsac1*Jsac1f';
sigbsacCAN1 = sqrt(varsac1f(1,1))
sigbsacUS1 = sqrt(varsac1f(2,2))
sigbsacav1 = sqrt(varsac1f(3,3))
tbsacCAN1 = abs(bsacCAN1 / sigbsacCAN1);
tbsacUS1 = abs(bsacUS1 / sigbsacUS1);
tbsacav1 = abs(bsacav1 / sigbsacav1);
pbsacCAN1 = 2*(1 - norm_cdf (tbsacCAN1))
pbsacUS1 = 2*(1 - norm_cdf (tbsacUS1))
pbsacav1 = 2*(1 - norm_cdf (tbsacav1))
%
%
xsac2 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij13 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressac2 = sacgrav(y,xsac2,W,W,info);
prt(ressac2,vnames);
[ressac2.beta
 ressac2.rho
 ressac2.lam]./ressac2.tstat
% Criteria information
logliksac2 = ressac2.lik;
AICsac2 = logliksac2*(2/n2) + ressac2.nvar*(2/n2)
SCsac2 = logliksac2*(2/n2) + ressac2.nvar*(log(n2)/n2)
% Borders
bsac2 = exp(inv(eye(n2,n2) - ressac2.rho*Wt)*(ressac2.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsacCANCAN2 = bsac2(1)
bsacCANUS2 = bsac2(11)
bsacUSUS2 = bsac2(411)
bsacUSCAN2 = bsac2(401)
bsacCAN2 = bsacCANCAN2/bsacCANUS2
bsacUS2 = bsacUSUS2/bsacUSCAN2
bsacav2 = sqrt(bsacCAN2*bsacUS2)
% p-value of partial border
Jsac2 = [(inv(eye(n2,n2) - ressac2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressac2.beta(4)*(inv(eye(n2,n2) - ressac2.rho*Wt))*Wt*(inv(eye(n2,n2) - ressac2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsac2 = [ressac2.minfo(4,4) ressac2.minfo(4,6)
    ressac2.minfo(6,4) ressac2.minfo(6,6)];
varsac2 = Jsac2*Vsac2*Jsac2';
sigbsacCANCAN2 = sqrt(varsac2(1,1))
sigbsacCANUS2 = sqrt(varsac2(11,11))
sigbsacUSUS2 = sqrt(varsac2(411,411))
sigbsacUSCAN2 = sqrt(varsac2(401,401))
tbsacCANCAN2 = abs(bsacCANCAN2 / sigbsacCANCAN2);
tbsacCANUS2 = abs(bsacCANUS2 / sigbsacCANUS2);
tbsacUSUS2 = abs(bsacUSUS2 / sigbsacUSUS2);
tbsacUSCAN2 = abs(bsacUSCAN2 / sigbsacUSCAN2);
pbsacCANCAN2 = 2*(1 - norm_cdf (tbsacCANCAN2))
pbsacCANUS2 = 2*(1 - norm_cdf (tbsacCANUS2))
pbsacUSUS2 = 2*(1 - norm_cdf (tbsacUSUS2))
pbsacUSCAN2 = 2*(1 - norm_cdf (tbsacUSCAN2))
% p-value of full border
Jsac2f = [1/bsacCANUS2  -(bsacCANCAN2)/((bsacCANUS2)^2)  0  0
    0  0  1/bsacUSCAN2  -(bsacUSUS2)/((bsacUSCAN2)^2)
    (1/(2*sqrt((bsacCANCAN2*bsacUSUS2)/(bsacCANUS2*bsacUSCAN2))))*((bsacUSUS2)/(bsacCANUS2*bsacUSCAN2))  -(1/(2*sqrt((bsacCANCAN2*bsacUSUS2)/(bsacCANUS2*bsacUSCAN2))))*((bsacCANCAN2*bsacUSUS2*bsacUSCAN2)/((bsacCANUS2*bsacUSCAN2)^2))  -(1/(2*sqrt((bsacCANCAN2*bsacUSUS2)/(bsacCANUS2*bsacUSCAN2))))*((bsacCANCAN2*bsacUSUS2*bsacCANUS2)/((bsacCANUS2*bsacUSCAN2)^2))  (1/(2*sqrt((bsacCANCAN2*bsacUSUS2)/(bsacCANUS2*bsacUSCAN2))))*((bsacCANCAN2)/(bsacCANUS2*bsacUSCAN2))];
subvarsac2 = [varsac2(1,1) varsac2(1,11) varsac2(1,411) varsac2(1,401)
              varsac2(11,1) varsac2(11,11) varsac2(11,411) varsac2(11,401)
              varsac2(411,1) varsac2(411,11) varsac2(411,411) varsac2(411,401)
              varsac2(401,1) varsac2(401,11) varsac2(401,411) varsac2(401,401)];
varsac2f = Jsac2f*subvarsac2*Jsac2f';
sigbsacCAN2 = sqrt(varsac2f(1,1))
sigbsacUS2 = sqrt(varsac2f(2,2))
sigbsacav2 = sqrt(varsac2f(3,3))
tbsacCAN2 = abs(bsacCAN2 / sigbsacCAN2);
tbsacUS2 = abs(bsacUS2 / sigbsacUS2);
tbsacav2 = abs(bsacav2 / sigbsacav2);
pbsacCAN2 = 2*(1 - norm_cdf (tbsacCAN2))
pbsacUS2 = 2*(1 - norm_cdf (tbsacUS2))
pbsacav2 = 2*(1 - norm_cdf (tbsacav2))
%
%
xsac3 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij23 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
vnames = strvcat('lnXij','constant','lndij','lnwage','VMborder','vm0');
ressac3 = sacgrav(y,xsac3,W,W,info);
prt(ressac3,vnames);
[ressac3.beta
 ressac3.rho
 ressac3.lam]./ressac3.tstat
% Criteria information
logliksac3 = ressac3.lik;
AICsac3 = logliksac3*(2/n2) + ressac3.nvar*(2/n2)
SCsac3 = logliksac3*(2/n2) + ressac3.nvar*(log(n2)/n2)
% Borders
bsac3 = exp(inv(eye(n2,n2) - ressac3.rho*Wt)*(ressac3.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsacCANCAN3 = bsac3(1)
bsacCANUS3 = bsac3(11)
bsacUSUS3 = bsac3(411)
bsacUSCAN3 = bsac3(401)
bsacCAN3 = bsacCANCAN3/bsacCANUS3
bsacUS3 = bsacUSUS3/bsacUSCAN3
bsacav3 = sqrt(bsacCAN3*bsacUS3)
% p-value of partial border
Jsac3 = [(inv(eye(n2,n2) - ressac3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressac3.beta(4)*(inv(eye(n2,n2) - ressac3.rho*Wt))*Wt*(inv(eye(n2,n2) - ressac3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsac3 = [ressac3.minfo(4,4) ressac3.minfo(4,6)
    ressac3.minfo(6,4) ressac3.minfo(6,6)];
varsac3 = Jsac3*Vsac3*Jsac3';
sigbsacCANCAN3 = sqrt(varsac3(1,1))
sigbsacCANUS3 = sqrt(varsac3(11,11))
sigbsacUSUS3 = sqrt(varsac3(411,411))
sigbsacUSCAN3 = sqrt(varsac3(401,401))
tbsacCANCAN3 = abs(bsacCANCAN3 / sigbsacCANCAN3);
tbsacCANUS3 = abs(bsacCANUS3 / sigbsacCANUS3);
tbsacUSUS3 = abs(bsacUSUS3 / sigbsacUSUS3);
tbsacUSCAN3 = abs(bsacUSCAN3 / sigbsacUSCAN3);
pbsacCANCAN3 = 2*(1 - norm_cdf (tbsacCANCAN3))
pbsacCANUS3 = 2*(1 - norm_cdf (tbsacCANUS3))
pbsacUSUS3 = 2*(1 - norm_cdf (tbsacUSUS3))
pbsacUSCAN3 = 2*(1 - norm_cdf (tbsacUSCAN3))
% p-value of full border
Jsac3f = [1/bsacCANUS3  -(bsacCANCAN3)/((bsacCANUS3)^2)  0  0
    0  0  1/bsacUSCAN3  -(bsacUSUS3)/((bsacUSCAN3)^2)
    (1/(2*sqrt((bsacCANCAN3*bsacUSUS3)/(bsacCANUS3*bsacUSCAN3))))*((bsacUSUS3)/(bsacCANUS3*bsacUSCAN3))  -(1/(2*sqrt((bsacCANCAN3*bsacUSUS3)/(bsacCANUS3*bsacUSCAN3))))*((bsacCANCAN3*bsacUSUS3*bsacUSCAN3)/((bsacCANUS3*bsacUSCAN3)^2))  -(1/(2*sqrt((bsacCANCAN3*bsacUSUS3)/(bsacCANUS3*bsacUSCAN3))))*((bsacCANCAN3*bsacUSUS3*bsacCANUS3)/((bsacCANUS3*bsacUSCAN3)^2))  (1/(2*sqrt((bsacCANCAN3*bsacUSUS3)/(bsacCANUS3*bsacUSCAN3))))*((bsacCANCAN3)/(bsacCANUS3*bsacUSCAN3))];
subvarsac3 = [varsac3(1,1) varsac3(1,11) varsac3(1,411) varsac3(1,401)
              varsac3(11,1) varsac3(11,11) varsac3(11,411) varsac3(11,401)
              varsac3(411,1) varsac3(411,11) varsac3(411,411) varsac3(411,401)
              varsac3(401,1) varsac3(401,11) varsac3(401,411) varsac3(401,401)];
varsac3f = Jsac3f*subvarsac3*Jsac3f';
sigbsacCAN3 = sqrt(varsac3f(1,1))
sigbsacUS3 = sqrt(varsac3f(2,2))
sigbsacav3 = sqrt(varsac3f(3,3))
tbsacCAN3 = abs(bsacCAN3 / sigbsacCAN3);
tbsacUS3 = abs(bsacUS3 / sigbsacUS3);
tbsacav3 = abs(bsacav3 / sigbsacav3);
pbsacCAN3 = 2*(1 - norm_cdf (tbsacCAN3))
pbsacUS3 = 2*(1 - norm_cdf (tbsacUS3))
pbsacav3 = 2*(1 - norm_cdf (tbsacav3))
%
% SARMA model
%
xsarma1 = [ones(n2,1) (eye(n2,n2) - Wt)*lndijAvW lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma1 = sarma(y,xsarma1,W);
% Criteria information
logliksarma1 = ressarma1.lik;
AICsarma1 = logliksarma1*(2/n2) + ressarma1.nvar*(2/n2)
SCsarma1 = logliksarma1*(2/n2) + ressarma1.nvar*(log(n2)/n2)
% Borders
bsarma1 = exp(inv(eye(n2,n2) - ressarma1.rho*Wt)*(ressarma1.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN1 = bsarma1(1)
bsarmaCANUS1 = bsarma1(11)
bsarmaUSUS1 = bsarma1(411)
bsarmaUSCAN1 = bsarma1(401)
bsarmaCAN1 = bsarmaCANCAN1/bsarmaCANUS1
bsarmaUS1 = bsarmaUSUS1/bsarmaUSCAN1
bsarmaav1 = sqrt(bsarmaCAN1*bsarmaUS1)
% p-value of partial border
Jsarma1 = [(inv(eye(n2,n2) - ressarma1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma1.beta(4)*(inv(eye(n2,n2) - ressarma1.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma1 = [ressarma1.minfo(4,4) ressarma1.minfo(4,6)
    ressarma1.minfo(6,4) ressarma1.minfo(6,6)];
varsarma1 = Jsarma1*Vsarma1*Jsarma1';
sigbsarmaCANCAN1 = sqrt(varsarma1(1,1))
sigbsarmaCANUS1 = sqrt(varsarma1(11,11))
sigbsarmaUSUS1 = sqrt(varsarma1(411,411))
sigbsarmaUSCAN1 = sqrt(varsarma1(401,401))
tbsarmaCANCAN1 = abs(bsarmaCANCAN1 / sigbsarmaCANCAN1);
tbsarmaCANUS1 = abs(bsarmaCANUS1 / sigbsarmaCANUS1);
tbsarmaUSUS1 = abs(bsarmaUSUS1 / sigbsarmaUSUS1);
tbsarmaUSCAN1 = abs(bsarmaUSCAN1 / sigbsarmaUSCAN1);
pbsarmaCANCAN1 = 2*(1 - norm_cdf (tbsarmaCANCAN1))
pbsarmaCANUS1 = 2*(1 - norm_cdf (tbsarmaCANUS1))
pbsarmaUSUS1 = 2*(1 - norm_cdf (tbsarmaUSUS1))
pbsarmaUSCAN1 = 2*(1 - norm_cdf (tbsarmaUSCAN1))
% p-value of full border
Jsarma1f = [1/bsarmaCANUS1  -(bsarmaCANCAN1)/((bsarmaCANUS1)^2)  0  0
    0  0  1/bsarmaUSCAN1  -(bsarmaUSUS1)/((bsarmaUSCAN1)^2)
    (1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))  -(1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1*bsarmaUSUS1*bsarmaUSCAN1)/((bsarmaCANUS1*bsarmaUSCAN1)^2))  -(1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1*bsarmaUSUS1*bsarmaCANUS1)/((bsarmaCANUS1*bsarmaUSCAN1)^2))  (1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1)/(bsarmaCANUS1*bsarmaUSCAN1))];
subvarsarma1 = [varsarma1(1,1) varsarma1(1,11) varsarma1(1,411) varsarma1(1,401)
              varsarma1(11,1) varsarma1(11,11) varsarma1(11,411) varsarma1(11,401)
              varsarma1(411,1) varsarma1(411,11) varsarma1(411,411) varsarma1(411,401)
              varsarma1(401,1) varsarma1(401,11) varsarma1(401,411) varsarma1(401,401)];
varsarma1f = Jsarma1f*subvarsarma1*Jsarma1f';
sigbsarmaCAN1 = sqrt(varsarma1f(1,1))
sigbsarmaUS1 = sqrt(varsarma1f(2,2))
sigbsarmaav1 = sqrt(varsarma1f(3,3))
tbsarmaCAN1 = abs(bsarmaCAN1 / sigbsarmaCAN1);
tbsarmaUS1 = abs(bsarmaUS1 / sigbsarmaUS1);
tbsarmaav1 = abs(bsarmaav1 / sigbsarmaav1);
pbsarmaCAN1 = 2*(1 - norm_cdf (tbsarmaCAN1))
pbsarmaUS1 = 2*(1 - norm_cdf (tbsarmaUS1))
pbsarmaav1 = 2*(1 - norm_cdf (tbsarmaav1))
%
%
xsarma2 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij13 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma2 = sarma(y,xsarma2,W);
% Criteria information
logliksarma2 = ressarma2.lik;
AICsarma2 = logliksarma2*(2/n2) + ressarma2.nvar*(2/n2)
SCsarma2 = logliksarma2*(2/n2) + ressarma2.nvar*(log(n2)/n2)
% Borders
bsarma2 = exp(inv(eye(n2,n2) - ressarma2.rho*Wt)*(ressarma2.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN2 = bsarma2(1)
bsarmaCANUS2 = bsarma2(11)
bsarmaUSUS2 = bsarma2(411)
bsarmaUSCAN2 = bsarma2(401)
bsarmaCAN2 = bsarmaCANCAN2/bsarmaCANUS2
bsarmaUS2 = bsarmaUSUS2/bsarmaUSCAN2
bsarmaav2 = sqrt(bsarmaCAN2*bsarmaUS2)
% p-value of partial border
Jsarma2 = [(inv(eye(n2,n2) - ressarma2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma2.beta(4)*(inv(eye(n2,n2) - ressarma2.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma2 = [ressarma2.minfo(4,4) ressarma2.minfo(4,6)
    ressarma2.minfo(6,4) ressarma2.minfo(6,6)];
varsarma2 = Jsarma2*Vsarma2*Jsarma2';
sigbsarmaCANCAN2 = sqrt(varsarma2(1,1))
sigbsarmaCANUS2 = sqrt(varsarma2(11,11))
sigbsarmaUSUS2 = sqrt(varsarma2(411,411))
sigbsarmaUSCAN2 = sqrt(varsarma2(401,401))
tbsarmaCANCAN2 = abs(bsarmaCANCAN2 / sigbsarmaCANCAN2);
tbsarmaCANUS2 = abs(bsarmaCANUS2 / sigbsarmaCANUS2);
tbsarmaUSUS2 = abs(bsarmaUSUS2 / sigbsarmaUSUS2);
tbsarmaUSCAN2 = abs(bsarmaUSCAN2 / sigbsarmaUSCAN2);
pbsarmaCANCAN2 = 2*(1 - norm_cdf (tbsarmaCANCAN2))
pbsarmaCANUS2 = 2*(1 - norm_cdf (tbsarmaCANUS2))
pbsarmaUSUS2 = 2*(1 - norm_cdf (tbsarmaUSUS2))
pbsarmaUSCAN2 = 2*(1 - norm_cdf (tbsarmaUSCAN2))
% p-value of full border
Jsarma2f = [1/bsarmaCANUS2  -(bsarmaCANCAN2)/((bsarmaCANUS2)^2)  0  0
    0  0  1/bsarmaUSCAN2  -(bsarmaUSUS2)/((bsarmaUSCAN2)^2)
    (1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))  -(1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2*bsarmaUSUS2*bsarmaUSCAN2)/((bsarmaCANUS2*bsarmaUSCAN2)^2))  -(1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2*bsarmaUSUS2*bsarmaCANUS2)/((bsarmaCANUS2*bsarmaUSCAN2)^2))  (1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2)/(bsarmaCANUS2*bsarmaUSCAN2))];
subvarsarma2 = [varsarma2(1,1) varsarma2(1,11) varsarma2(1,411) varsarma2(1,401)
              varsarma2(11,1) varsarma2(11,11) varsarma2(11,411) varsarma2(11,401)
              varsarma2(411,1) varsarma2(411,11) varsarma2(411,411) varsarma2(411,401)
              varsarma2(401,1) varsarma2(401,11) varsarma2(401,411) varsarma2(401,401)];
varsarma2f = Jsarma2f*subvarsarma2*Jsarma2f';
sigbsarmaCAN2 = sqrt(varsarma2f(1,1))
sigbsarmaUS2 = sqrt(varsarma2f(2,2))
sigbsarmaav2 = sqrt(varsarma2f(3,3))
tbsarmaCAN2 = abs(bsarmaCAN2 / sigbsarmaCAN2);
tbsarmaUS2 = abs(bsarmaUS2 / sigbsarmaUS2);
tbsarmaav2 = abs(bsarmaav2 / sigbsarmaav2);
pbsarmaCAN2 = 2*(1 - norm_cdf (tbsarmaCAN2))
pbsarmaUS2 = 2*(1 - norm_cdf (tbsarmaUS2))
pbsarmaav2 = 2*(1 - norm_cdf (tbsarmaav2))
%
%
xsarma3 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij23 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma3 = sarma(y,xsarma3,W);
% Criteria information
logliksarma3 = ressarma3.lik;
AICsarma3 = logliksarma3*(2/n2) + ressarma3.nvar*(2/n2)
SCsarma3 = logliksarma3*(2/n2) + ressarma3.nvar*(log(n2)/n2)
% Borders
bsarma3 = exp(inv(eye(n2,n2) - ressarma3.rho*Wt)*(ressarma3.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN3 = bsarma3(1)
bsarmaCANUS3 = bsarma3(11)
bsarmaUSUS3 = bsarma3(411)
bsarmaUSCAN3 = bsarma3(401)
bsarmaCAN3 = bsarmaCANCAN3/bsarmaCANUS3
bsarmaUS3 = bsarmaUSUS3/bsarmaUSCAN3
bsarmaav3 = sqrt(bsarmaCAN3*bsarmaUS3)
% p-value of partial border
Jsarma3 = [(inv(eye(n2,n2) - ressarma3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma3.beta(4)*(inv(eye(n2,n2) - ressarma3.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma3 = [ressarma3.minfo(4,4) ressarma3.minfo(4,6)
    ressarma3.minfo(6,4) ressarma3.minfo(6,6)];
varsarma3 = Jsarma3*Vsarma3*Jsarma3';
sigbsarmaCANCAN3 = sqrt(varsarma3(1,1))
sigbsarmaCANUS3 = sqrt(varsarma3(11,11))
sigbsarmaUSUS3 = sqrt(varsarma3(411,411))
sigbsarmaUSCAN3 = sqrt(varsarma3(401,401))
tbsarmaCANCAN3 = abs(bsarmaCANCAN3 / sigbsarmaCANCAN3);
tbsarmaCANUS3 = abs(bsarmaCANUS3 / sigbsarmaCANUS3);
tbsarmaUSUS3 = abs(bsarmaUSUS3 / sigbsarmaUSUS3);
tbsarmaUSCAN3 = abs(bsarmaUSCAN3 / sigbsarmaUSCAN3);
pbsarmaCANCAN3 = 2*(1 - norm_cdf (tbsarmaCANCAN3))
pbsarmaCANUS3 = 2*(1 - norm_cdf (tbsarmaCANUS3))
pbsarmaUSUS3 = 2*(1 - norm_cdf (tbsarmaUSUS3))
pbsarmaUSCAN3 = 2*(1 - norm_cdf (tbsarmaUSCAN3))
% p-value of full border
Jsarma3f = [1/bsarmaCANUS3  -(bsarmaCANCAN3)/((bsarmaCANUS3)^2)  0  0
    0  0  1/bsarmaUSCAN3  -(bsarmaUSUS3)/((bsarmaUSCAN3)^2)
    (1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))  -(1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3*bsarmaUSUS3*bsarmaUSCAN3)/((bsarmaCANUS3*bsarmaUSCAN3)^2))  -(1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3*bsarmaUSUS3*bsarmaCANUS3)/((bsarmaCANUS3*bsarmaUSCAN3)^2))  (1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3)/(bsarmaCANUS3*bsarmaUSCAN3))];
subvarsarma3 = [varsarma3(1,1) varsarma3(1,11) varsarma3(1,411) varsarma3(1,401)
              varsarma3(11,1) varsarma3(11,11) varsarma3(11,411) varsarma3(11,401)
              varsarma3(411,1) varsarma3(411,11) varsarma3(411,411) varsarma3(411,401)
              varsarma3(401,1) varsarma3(401,11) varsarma3(401,411) varsarma3(401,401)];
varsarma3f = Jsarma3f*subvarsarma3*Jsarma3f';
sigbsarmaCAN3 = sqrt(varsarma3f(1,1))
sigbsarmaUS3 = sqrt(varsarma3f(2,2))
sigbsarmaav3 = sqrt(varsarma3f(3,3))
tbsarmaCAN3 = abs(bsarmaCAN3 / sigbsarmaCAN3);
tbsarmaUS3 = abs(bsarmaUS3 / sigbsarmaUS3);
tbsarmaav3 = abs(bsarmaav3 / sigbsarmaav3);
pbsarmaCAN3 = 2*(1 - norm_cdf (tbsarmaCAN3))
pbsarmaUS3 = 2*(1 - norm_cdf (tbsarmaUS3))
pbsarmaav3 = 2*(1 - norm_cdf (tbsarmaav3))
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Robustness (US-US flow*1.3)
y13 = lnXij13 - lnYi - lnYj;
%
% SARMA model
%
xsarma1 = [ones(n2,1) (eye(n2,n2) - Wt)*lndijAvW lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma1 = sarma(y13,xsarma1,W);
% Criteria information
logliksarma1 = ressarma1.lik;
AICsarma1 = logliksarma1*(2/n2) + ressarma1.nvar*(2/n2)
SCsarma1 = logliksarma1*(2/n2) + ressarma1.nvar*(log(n2)/n2)
% Borders
bsarma1 = exp(inv(eye(n2,n2) - ressarma1.rho*Wt)*(ressarma1.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN1 = bsarma1(1)
bsarmaCANUS1 = bsarma1(11)
bsarmaUSUS1 = bsarma1(411)
bsarmaUSCAN1 = bsarma1(401)
bsarmaCAN1 = bsarmaCANCAN1/bsarmaCANUS1
bsarmaUS1 = bsarmaUSUS1/bsarmaUSCAN1
bsarmaav1 = sqrt(bsarmaCAN1*bsarmaUS1)
% p-value of partial border
Jsarma1 = [(inv(eye(n2,n2) - ressarma1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma1.beta(4)*(inv(eye(n2,n2) - ressarma1.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma1.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma1 = [ressarma1.minfo(4,4) ressarma1.minfo(4,6)
    ressarma1.minfo(6,4) ressarma1.minfo(6,6)];
varsarma1 = Jsarma1*Vsarma1*Jsarma1';
sigbsarmaCANCAN1 = sqrt(varsarma1(1,1))
sigbsarmaCANUS1 = sqrt(varsarma1(11,11))
sigbsarmaUSUS1 = sqrt(varsarma1(411,411))
sigbsarmaUSCAN1 = sqrt(varsarma1(401,401))
tbsarmaCANCAN1 = abs(bsarmaCANCAN1 / sigbsarmaCANCAN1);
tbsarmaCANUS1 = abs(bsarmaCANUS1 / sigbsarmaCANUS1);
tbsarmaUSUS1 = abs(bsarmaUSUS1 / sigbsarmaUSUS1);
tbsarmaUSCAN1 = abs(bsarmaUSCAN1 / sigbsarmaUSCAN1);
pbsarmaCANCAN1 = 2*(1 - norm_cdf (tbsarmaCANCAN1))
pbsarmaCANUS1 = 2*(1 - norm_cdf (tbsarmaCANUS1))
pbsarmaUSUS1 = 2*(1 - norm_cdf (tbsarmaUSUS1))
pbsarmaUSCAN1 = 2*(1 - norm_cdf (tbsarmaUSCAN1))
% p-value of full border
Jsarma1f = [1/bsarmaCANUS1  -(bsarmaCANCAN1)/((bsarmaCANUS1)^2)  0  0
    0  0  1/bsarmaUSCAN1  -(bsarmaUSUS1)/((bsarmaUSCAN1)^2)
    (1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))  -(1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1*bsarmaUSUS1*bsarmaUSCAN1)/((bsarmaCANUS1*bsarmaUSCAN1)^2))  -(1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1*bsarmaUSUS1*bsarmaCANUS1)/((bsarmaCANUS1*bsarmaUSCAN1)^2))  (1/(2*sqrt((bsarmaCANCAN1*bsarmaUSUS1)/(bsarmaCANUS1*bsarmaUSCAN1))))*((bsarmaCANCAN1)/(bsarmaCANUS1*bsarmaUSCAN1))];
subvarsarma1 = [varsarma1(1,1) varsarma1(1,11) varsarma1(1,411) varsarma1(1,401)
              varsarma1(11,1) varsarma1(11,11) varsarma1(11,411) varsarma1(11,401)
              varsarma1(411,1) varsarma1(411,11) varsarma1(411,411) varsarma1(411,401)
              varsarma1(401,1) varsarma1(401,11) varsarma1(401,411) varsarma1(401,401)];
varsarma1f = Jsarma1f*subvarsarma1*Jsarma1f';
sigbsarmaCAN1 = sqrt(varsarma1f(1,1))
sigbsarmaUS1 = sqrt(varsarma1f(2,2))
sigbsarmaav1 = sqrt(varsarma1f(3,3))
tbsarmaCAN1 = abs(bsarmaCAN1 / sigbsarmaCAN1);
tbsarmaUS1 = abs(bsarmaUS1 / sigbsarmaUS1);
tbsarmaav1 = abs(bsarmaav1 / sigbsarmaav1);
pbsarmaCAN1 = 2*(1 - norm_cdf (tbsarmaCAN1))
pbsarmaUS1 = 2*(1 - norm_cdf (tbsarmaUS1))
pbsarmaav1 = 2*(1 - norm_cdf (tbsarmaav1))
%
%
xsarma2 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij13 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma2 = sarma(y13,xsarma2,W);
% Criteria information
logliksarma2 = ressarma2.lik;
AICsarma2 = logliksarma2*(2/n2) + ressarma2.nvar*(2/n2)
SCsarma2 = logliksarma2*(2/n2) + ressarma2.nvar*(log(n2)/n2)
% Borders
bsarma2 = exp(inv(eye(n2,n2) - ressarma2.rho*Wt)*(ressarma2.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN2 = bsarma2(1)
bsarmaCANUS2 = bsarma2(11)
bsarmaUSUS2 = bsarma2(411)
bsarmaUSCAN2 = bsarma2(401)
bsarmaCAN2 = bsarmaCANCAN2/bsarmaCANUS2
bsarmaUS2 = bsarmaUSUS2/bsarmaUSCAN2
bsarmaav2 = sqrt(bsarmaCAN2*bsarmaUS2)
% p-value of partial border
Jsarma2 = [(inv(eye(n2,n2) - ressarma2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma2.beta(4)*(inv(eye(n2,n2) - ressarma2.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma2.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma2 = [ressarma2.minfo(4,4) ressarma2.minfo(4,6)
    ressarma2.minfo(6,4) ressarma2.minfo(6,6)];
varsarma2 = Jsarma2*Vsarma2*Jsarma2';
sigbsarmaCANCAN2 = sqrt(varsarma2(1,1))
sigbsarmaCANUS2 = sqrt(varsarma2(11,11))
sigbsarmaUSUS2 = sqrt(varsarma2(411,411))
sigbsarmaUSCAN2 = sqrt(varsarma2(401,401))
tbsarmaCANCAN2 = abs(bsarmaCANCAN2 / sigbsarmaCANCAN2);
tbsarmaCANUS2 = abs(bsarmaCANUS2 / sigbsarmaCANUS2);
tbsarmaUSUS2 = abs(bsarmaUSUS2 / sigbsarmaUSUS2);
tbsarmaUSCAN2 = abs(bsarmaUSCAN2 / sigbsarmaUSCAN2);
pbsarmaCANCAN2 = 2*(1 - norm_cdf (tbsarmaCANCAN2))
pbsarmaCANUS2 = 2*(1 - norm_cdf (tbsarmaCANUS2))
pbsarmaUSUS2 = 2*(1 - norm_cdf (tbsarmaUSUS2))
pbsarmaUSCAN2 = 2*(1 - norm_cdf (tbsarmaUSCAN2))
% p-value of full border
Jsarma2f = [1/bsarmaCANUS2  -(bsarmaCANCAN2)/((bsarmaCANUS2)^2)  0  0
    0  0  1/bsarmaUSCAN2  -(bsarmaUSUS2)/((bsarmaUSCAN2)^2)
    (1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))  -(1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2*bsarmaUSUS2*bsarmaUSCAN2)/((bsarmaCANUS2*bsarmaUSCAN2)^2))  -(1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2*bsarmaUSUS2*bsarmaCANUS2)/((bsarmaCANUS2*bsarmaUSCAN2)^2))  (1/(2*sqrt((bsarmaCANCAN2*bsarmaUSUS2)/(bsarmaCANUS2*bsarmaUSCAN2))))*((bsarmaCANCAN2)/(bsarmaCANUS2*bsarmaUSCAN2))];
subvarsarma2 = [varsarma2(1,1) varsarma2(1,11) varsarma2(1,411) varsarma2(1,401)
              varsarma2(11,1) varsarma2(11,11) varsarma2(11,411) varsarma2(11,401)
              varsarma2(411,1) varsarma2(411,11) varsarma2(411,411) varsarma2(411,401)
              varsarma2(401,1) varsarma2(401,11) varsarma2(401,411) varsarma2(401,401)];
varsarma2f = Jsarma2f*subvarsarma2*Jsarma2f';
sigbsarmaCAN2 = sqrt(varsarma2f(1,1))
sigbsarmaUS2 = sqrt(varsarma2f(2,2))
sigbsarmaav2 = sqrt(varsarma2f(3,3))
tbsarmaCAN2 = abs(bsarmaCAN2 / sigbsarmaCAN2);
tbsarmaUS2 = abs(bsarmaUS2 / sigbsarmaUS2);
tbsarmaav2 = abs(bsarmaav2 / sigbsarmaav2);
pbsarmaCAN2 = 2*(1 - norm_cdf (tbsarmaCAN2))
pbsarmaUS2 = 2*(1 - norm_cdf (tbsarmaUS2))
pbsarmaav2 = 2*(1 - norm_cdf (tbsarmaav2))
%
%
xsarma3 = [ones(n2,1) (eye(n2,n2) - Wt)*lndij23 lnwage (eye(n2,n2) - Wt)*(vmCAN+vmUS) vm0];
ressarma3 = sarma(y13,xsarma3,W);
% Criteria information
logliksarma3 = ressarma3.lik;
AICsarma3 = logliksarma3*(2/n2) + ressarma3.nvar*(2/n2)
SCsarma3 = logliksarma3*(2/n2) + ressarma3.nvar*(log(n2)/n2)
% Borders
bsarma3 = exp(inv(eye(n2,n2) - ressarma3.rho*Wt)*(ressarma3.beta(4)*(eye(n2,n2) - Wt)*(vmCAN+vmUS)));
bsarmaCANCAN3 = bsarma3(1)
bsarmaCANUS3 = bsarma3(11)
bsarmaUSUS3 = bsarma3(411)
bsarmaUSCAN3 = bsarma3(401)
bsarmaCAN3 = bsarmaCANCAN3/bsarmaCANUS3
bsarmaUS3 = bsarmaUSUS3/bsarmaUSCAN3
bsarmaav3 = sqrt(bsarmaCAN3*bsarmaUS3)
% p-value of partial border
Jsarma3 = [(inv(eye(n2,n2) - ressarma3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS) ressarma3.beta(4)*(inv(eye(n2,n2) - ressarma3.rho*Wt))*Wt*(inv(eye(n2,n2) - ressarma3.rho*Wt))*(eye(n2,n2) - Wt)*(vmCAN+vmUS)];
Vsarma3 = [ressarma3.minfo(4,4) ressarma3.minfo(4,6)
    ressarma3.minfo(6,4) ressarma3.minfo(6,6)];
varsarma3 = Jsarma3*Vsarma3*Jsarma3';
sigbsarmaCANCAN3 = sqrt(varsarma3(1,1))
sigbsarmaCANUS3 = sqrt(varsarma3(11,11))
sigbsarmaUSUS3 = sqrt(varsarma3(411,411))
sigbsarmaUSCAN3 = sqrt(varsarma3(401,401))
tbsarmaCANCAN3 = abs(bsarmaCANCAN3 / sigbsarmaCANCAN3);
tbsarmaCANUS3 = abs(bsarmaCANUS3 / sigbsarmaCANUS3);
tbsarmaUSUS3 = abs(bsarmaUSUS3 / sigbsarmaUSUS3);
tbsarmaUSCAN3 = abs(bsarmaUSCAN3 / sigbsarmaUSCAN3);
pbsarmaCANCAN3 = 2*(1 - norm_cdf (tbsarmaCANCAN3))
pbsarmaCANUS3 = 2*(1 - norm_cdf (tbsarmaCANUS3))
pbsarmaUSUS3 = 2*(1 - norm_cdf (tbsarmaUSUS3))
pbsarmaUSCAN3 = 2*(1 - norm_cdf (tbsarmaUSCAN3))
% p-value of full border
Jsarma3f = [1/bsarmaCANUS3  -(bsarmaCANCAN3)/((bsarmaCANUS3)^2)  0  0
    0  0  1/bsarmaUSCAN3  -(bsarmaUSUS3)/((bsarmaUSCAN3)^2)
    (1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))  -(1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3*bsarmaUSUS3*bsarmaUSCAN3)/((bsarmaCANUS3*bsarmaUSCAN3)^2))  -(1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3*bsarmaUSUS3*bsarmaCANUS3)/((bsarmaCANUS3*bsarmaUSCAN3)^2))  (1/(2*sqrt((bsarmaCANCAN3*bsarmaUSUS3)/(bsarmaCANUS3*bsarmaUSCAN3))))*((bsarmaCANCAN3)/(bsarmaCANUS3*bsarmaUSCAN3))];
subvarsarma3 = [varsarma3(1,1) varsarma3(1,11) varsarma3(1,411) varsarma3(1,401)
              varsarma3(11,1) varsarma3(11,11) varsarma3(11,411) varsarma3(11,401)
              varsarma3(411,1) varsarma3(411,11) varsarma3(411,411) varsarma3(411,401)
              varsarma3(401,1) varsarma3(401,11) varsarma3(401,411) varsarma3(401,401)];
varsarma3f = Jsarma3f*subvarsarma3*Jsarma3f';
sigbsarmaCAN3 = sqrt(varsarma3f(1,1))
sigbsarmaUS3 = sqrt(varsarma3f(2,2))
sigbsarmaav3 = sqrt(varsarma3f(3,3))
tbsarmaCAN3 = abs(bsarmaCAN3 / sigbsarmaCAN3);
tbsarmaUS3 = abs(bsarmaUS3 / sigbsarmaUS3);
tbsarmaav3 = abs(bsarmaav3 / sigbsarmaav3);
pbsarmaCAN3 = 2*(1 - norm_cdf (tbsarmaCAN3))
pbsarmaUS3 = 2*(1 - norm_cdf (tbsarmaUS3))
pbsarmaav3 = 2*(1 - norm_cdf (tbsarmaav3))
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SARMA beta hetero + 2 rho + 2 lambda
%
dcan = data(:,16);
dus = data(:,17);
%
%
% Variables
yc = lnXij - lnYi -lnYj + log(Sp);
%
zh = [lnwage (eye(n2,n2) - Wt)*lndij13 (eye(n2,n2) - Wt)*(vmCAN+vmUS)];
%
V = zeros(n2,n);
for j=1:n;
    V((j-1)*n+1:j*n,j) = 1;
end;
%
Z1sarma = zeros(n2,3*n);
for i=0:n-1;
Z1sarma(:,((i*3)+1):((i+1)*3)) = diag(V(:,i+1))*zh;
end;
%
Zsarma = [Z1sarma vm0];
%
resh = sarmadouble(yc,Zsarma,W,dcan,dus)
%
constanti = zeros(40,1);
for i=0:39
constanti(i+1) = resh.beta((i+1)*3-2);
end
disti = zeros(40,1);
for i=0:39
disti(i+1) = resh.beta((i+1)*3-1);
end
thetai = zeros(40,1);
for i=0:39
thetai(i+1) = resh.beta((i+1)*3);
end
pconstanti = zeros(40,1);
for i=0:39
pconstanti(i+1) = resh.pvalue((i+1)*3-2);
end
pdisti = zeros(40,1);
for i=0:39
pdisti(i+1) = resh.pvalue((i+1)*3-1);
end
pthetai = zeros(40,1);
for i=0:39
pthetai(i+1) = resh.pvalue((i+1)*3);
end
%
rho = [resh.rhocan(1,1)*ones(400,1)
    resh.rhous(1,1)*ones(1200,1)];
%
Wcan = sparse(diag(dcan)*W);
Wus = sparse(diag(dus)*W);
A = inv(eye(n2,n2) - resh.rhocan(1,1)*Wcan - resh.rhous(1,1)*Wus);
%
b = (eye(n2,n2) - Wt)*(vmCAN+vmUS);
%
% Border
%
th = zeros(1600,1);
for i=0:39
   th(i*40+1:40*(i+1)) = thetai(i+1); 
end
%
border = exp(A*(th.*b));
%
% Partial Borders
bsarmahCANCAN = zeros(10,1);
for i=0:9
bsarmahCANCAN(i+1) = border((i+1)*40-39);
end
bsarmahCANUS = zeros(10,1);
for i=0:9
bsarmahCANUS(i+1) = border((i+1)*40-29);
end
bsarmahUSUS = zeros(30,1);
for i=10:39
bsarmahUSUS(i-9) = border((i+1)*40-29);
end
bsarmahUSCAN = zeros(30,1);
for i=10:39
bsarmahUSCAN(i-9) = border((i+1)*40-39);
end
borderintra = [bsarmahCANCAN
               bsarmahUSUS];
borderinter = [bsarmahCANUS
               bsarmahUSCAN];
%
% Full Borders
bsarmahCAN = bsarmahCANCAN./bsarmahCANUS;
bsarmahUS = bsarmahUSUS./bsarmahUSCAN;
fullborder = [bsarmahCAN
              bsarmahUS];
%
% Local impact of D
%
dist = zeros(1600,1);
for i=0:39
   dist(i*40+1:40*(i+1)) = disti(i+1); 
end
%
distance = diag(diag(dist)*(A*(eye(n2,n2) - Wt)));
%
distloc = zeros(40,1);
for i=0:39
distloc(i+1) = distance(i*40+1);
end
%
% Results
%
results = [constanti pconstanti disti pdisti thetai pthetai distloc borderintra borderinter fullborder];
%
%
