%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
load data91.txt;
xlong = data91(:,1);
ylat = data91(:,2);
lny60 = data91(:,3);
lny95 = data91(:,4);
gy = data91(:,5);
lns = data91(:,6);
lnngd = data91(:,7);
n = length(data91);

info.lflag = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial weight matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% radian
X = (pi/180)*ylat;
Y = (pi/180)*xlong;

% Distance arc matrix (D)
T1=repmat(X',n,1);
T2=repmat(Y',n,1);
T3=repmat(X,1,n);
T4=repmat(Y,1,n);
D = (6378 * acos(cos(abs(T2-T4)).*cos(T1).*cos(T3) + sin(T1).*sin(T3))).*(ones(n,n)-eye(n,n));

% Inverse square distance matrix 
Wd2 = normw((ones(n)./(D+eye(n))).^2 - eye(n));

% Inverse square exponential matrix
We2 = normw(exp(-2*(D/1000)) - eye(n));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unconditional convergence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xuc = [ones(n,1) lny60];
vuc = strvcat('gy','constant','lny60');
resuc = ols(gy,xuc);
prt(resuc,vuc);

% Structural parameters

constant = resuc.beta(1);
a = resuc.beta(2);
siguc = resuc.sige;
t = 35;
lambdauc = - (1/t) * log(1+a*t)

% Spatial autocorrelation tests

moranucd2 = moran(gy,xuc,Wd2);
prt(moranucd2);

moranuce2 = moran(gy,xuc,We2);
prt(moranuce2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solow model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%% Level Solow unconstrained

xsolow = [ones(n,1) lns lnngd];
vsolow = strvcat('lny95','constant','lns','lnngd');
ressolow = ols(lny95,xsolow);
prt(ressolow,vsolow);

% Moran I test

moransolowd2 = moran(lny95,xsolow,Wd2);
prt(moransolowd2);

moransolowe2 = moran(lny95,xsolow,We2);
prt(moransolowe2);

%%%%%%%%%%%%%%%%%%%%%% Level Solow constrained

xsolowc = [ones(n,1) lns-lnngd];
vsolowc = strvcat('lny95','constant','lnsngd');
ressolowc = ols(lny95,xsolowc);
prt(ressolowc,vsolowc);

% Structural parameters

constant = ressolowc.beta(1);
a = ressolowc.beta(2);
alphasolowc = a / (1 + a)
sigsolowc = ressolowc.sige;
Jsolowc = 1 / (1 + a)^2;
sigalphasolowc = sigsolowc * Jsolowc;
talphasolowc = alphasolowc / sigalphasolowc;
palphasolowc = 2*(1 - norm_cdf (talphasolowc))

% Wald test of theoretical contraints

waldf(ressolowc,ressolow);

% Moran I test

moransolowcd2 = moran(lny95,xsolowc,Wd2);
prt(moransolowcd2);

moransolowce2 = moran(lny95,xsolowc,We2);
prt(moransolowce2);

%%%%%%%%%%%%%%%%%%%%%% Convergence Solow

xsolowconv = [ones(n,1) lny60 lns lnngd];
vsolowconv = strvcat('gy','constant','lny60','lns','lnngd');
ressolowconv = ols(gy,xsolowconv);
prt(ressolowconv,vsolowconv);

% Structural parameters

constant = ressolowconv.beta(1);
a = ressolowconv.beta(2);
sigsolowconv = ressolowconv.sige;
t = 35;
lambdasolowconv = - (1/t) * log(1+a*t)
Jsolowconv = - 1 / (1+a*t);
siglambdasolowconv = abs(sigsolowconv*Jsolowconv);
tlambdasolowconv = abs(lambdasolowconv / siglambdasolowconv);
plambdasolowconv = 2*(1 - norm_cdf (tlambdasolowconv))

% Moran I test

moransolowconvd2 = moran(gy,xsolowconv,Wd2);
prt(moransolowconvd2);

moransolowconve2 = moran(gy,xsolowconv,We2);
prt(moransolowconve2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial Solow SEM (spatial error model)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%% Level Solow SEM unconstrained

xsolow = [ones(n,1) lns lnngd];
vsolow = strvcat('lny95','constant','lns','lnngd');

semsolowd2 = sem(lny95,xsolow,Wd2,info);
prt(semsolowd2,vsolow);
semsolowe2 = sem(lny95,xsolow,We2,info);
prt(semsolowe2,vsolow);

%%%%%%%%%%%%%%%%%%%%%% Level Solow SEM constrained

xsolowc = [ones(n,1) lns-lnngd];
vsolowc = strvcat('lny95','constant','lnsngd');

semsolowcd2 = sem(lny95,xsolowc,Wd2,info);
prt(semsolowcd2,vsolowc);
semsolowce2 = sem(lny95,xsolowc,We2,info);
prt(semsolowce2,vsolowc);

% Structural parameters (Wd2)

constant = semsolowcd2.beta(1);
a = semsolowcd2.beta(2);
alphasemsolowcd2 = a / (1 + a)
sigsemsolowcd2 = semsolowcd2.sige;
Jsemsolowcd2 = 1 / (1 + a)^2;
sigalphasemsolowcd2 = sigsemsolowcd2 * Jsemsolowcd2;
talphasemsolowcd2 = alphasemsolowcd2 / sigalphasemsolowcd2;
palphasemsolowcd2 = 2*(1 - norm_cdf (talphasemsolowcd2))

% Structural parameters (We2)

constant = semsolowce2.beta(1);
a = semsolowce2.beta(2);
alphasemsolowce2 = a / (1 + a)
sigsemsolowce2 = semsolowce2.sige;
Jsemsolowce2 = 1 / (1 + a)^2;
sigalphasemsolowce2 = sigsemsolowce2 * Jsemsolowce2;
talphasemsolowce2 = alphasemsolowce2 / sigalphasemsolowce2;
palphasemsolowce2 = 2*(1 - norm_cdf (talphasemsolowce2))

% LR test for theoretical constraint

logliksemsolowd2 = semsolowd2.lik;
logliksemsolowcd2 = semsolowcd2.lik;
LRsemsolowd2 = -2*(logliksemsolowcd2 - logliksemsolowd2)
psemsolowd2 = 1-chis_prb(LRsemsolowd2,1)

logliksemsolowe2 = semsolowe2.lik;
logliksemsolowce2 = semsolowce2.lik;
LRsemsolowe2 = -2*(logliksemsolowce2 - logliksemsolowe2)
psemsolowe2 = 1-chis_prb(LRsemsolowe2,1)

%%%%%%%%%%%%%%%%%%%%%% Convergence Solow SEM

xsolowconv = [ones(n,1) lny60 lns lnngd];
vsolowconv = strvcat('gy','constant','lny60','lns','lnngd');

semsolowconvd2 = sem(gy,xsolowconv,Wd2,info);
prt(semsolowconvd2,vsolowconv);
semsolowconve2 = sem(gy,xsolowconv,We2,info);
prt(semsolowconve2,vsolowconv);

% Structural parameters (Wd2)

constant = semsolowconvd2.beta(1);
a = semsolowconvd2.beta(2);
sigsemsolowconvd2 = semsolowconvd2.sige;
t = 35;
lambdasemsolowconvd2 = - (1/t) * log(1+a*t)
Jsemsolowconvd2 = - 1 / (1+a*t);
sigsemsolowconvd2 = abs(sigsemsolowconvd2*Jsemsolowconvd2);
tlambdasemsolowconvd2 = abs(lambdasemsolowconvd2 / sigsemsolowconvd2);
plambdasemsolowconvd2 = 2*(1 - norm_cdf (tlambdasemsolowconvd2))

% Structural parameters (We2)

constant = semsolowconve2.beta(1);
a = semsolowconve2.beta(2);
sigsemsolowconve2 = semsolowconve2.sige;
t = 35;
lambdasemsolowconve2 = - (1/t) * log(1+a*t)
Jsemsolowconve2 = - 1 / (1+a*t);
sigsemsolowconve2 = abs(sigsemsolowconve2*Jsemsolowconve2);
tlambdasemsolowconve2 = abs(lambdasemsolowconve2 / sigsemsolowconve2);
plambdasemsolowconve2 = 2*(1 - norm_cdf (tlambdasemsolowconve2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial Solow SDM (spatial Durbin model)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%% Level Solow SDM unconstrained

xsolow = [ones(n,1) lns lnngd];
vsolow = strvcat('lny95','constant','lns','lnngd');

sdmsolowd2 = sdm(lny95,xsolow,Wd2,info);
prt(sdmsolowd2,vsolow);
sdmsolowe2 = sdm(lny95,xsolow,We2,info);
prt(sdmsolowe2,vsolow);

% LR test for common factor

logliksemsolowd2 = semsolowd2.lik;
logliksdmsolowd2 = sdmsolowd2.lik;
LRsolowd2 = -2*(logliksemsolowd2 - logliksdmsolowd2)
psolowd2 = 1-chis_prb(LRsolowd2,2)

logliksemsolowe2 = semsolowe2.lik;
logliksdmsolowe2 = sdmsolowe2.lik;
LRsolowe2 = -2*(logliksemsolowe2 - logliksdmsolowe2)
psolowe2 = 1-chis_prb(LRsolowe2,2)

% Elasticities, lns (Wd2)

constant = sdmsolowd2.beta(1);
a = sdmsolowd2.beta(2);
b = sdmsolowd2.beta(4);
c = sdmsolowd2.rho;
sigsdmsolowd2 = sdmsolowd2.sige;

Esdmsolowd2lns = (inv(eye(n,n) - c*Wd2)) * (a*eye(n,n) + b*Wd2);

% Elasticities, lnngd (Wd2)

constant = sdmsolowd2.beta(1);
a = sdmsolowd2.beta(3);
b = sdmsolowd2.beta(5);
c = sdmsolowd2.rho;
sigsdmsolowd2 = sdmsolowd2.sige;

Esdmsolowd2lnngd = (inv(eye(n,n) - c*Wd2)) * (a*eye(n,n) + b*Wd2);

% Elasticities, lns (We2)

constant = sdmsolowe2.beta(1);
a = sdmsolowe2.beta(2);
b = sdmsolowe2.beta(4);
c = sdmsolowe2.rho;
sigsdmsolowe2 = sdmsolowe2.sige;

Esdmsolowe2lns = (inv(eye(n,n) - c*We2)) * (a*eye(n,n) + b*We2);

% Elasticities, lnngd (We2)

constant = sdmsolowe2.beta(1);
a = sdmsolowe2.beta(3);
b = sdmsolowe2.beta(5);
c = sdmsolowe2.rho;
sigsdmsolowe2 = sdmsolowe2.sige;

Esdmsolowe2lnngd = (inv(eye(n,n) - c*We2)) * (a*eye(n,n) + b*We2);

%%%%%%%%%%%%%%%%%%%%%% Level Solow SDM constrained

xsolowc = [ones(n,1) lns-lnngd];
vsolowc = strvcat('lny95','constant','lnsngd');

sdmsolowcd2 = sdm(lny95,xsolowc,Wd2,info);
prt(sdmsolowcd2,vsolowc);
sdmsolowce2 = sdm(lny95,xsolowc,We2,info);
prt(sdmsolowce2,vsolowc);

% Structural parameters (Wd2)

constant = sdmsolowcd2.beta(1);
a = sdmsolowcd2.beta(2);
b = sdmsolowcd2.beta(3);
c = sdmsolowcd2.rho;

alphasdmsolowcd2 = b / (b - c)
phisdmsolowcd2 = a / (1 + a) - alphasdmsolowcd2
gammasdmsolowcd2 = - (b - c) / (1 + a)
REsdmsolowcd2 = (a + b) / (1 + a + b - c)
Jsdmsolowcd2=zeros(4,3);
Jsdmsolowcd2(1,1) = 0;
Jsdmsolowcd2(1,2) = -c / (b - c)^2;
Jsdmsolowcd2(1,3) = b / (b - c)^2;
Jsdmsolowcd2(2,1) = 1 / (1 + a)^2;
Jsdmsolowcd2(2,2) = c / (b - c)^2;
Jsdmsolowcd2(2,3) = - b / (b - c)^2;
Jsdmsolowcd2(3,1) = (b - c) / (1 + a)^2;
Jsdmsolowcd2(3,2) = - 1 / (1 + a);
Jsdmsolowcd2(3,3) = 1 / (1 + a);
Jsdmsolowcd2(4,1) = (1 - c) / (1 + a + b - c)^2;
Jsdmsolowcd2(4,2) = (1 - c) / (1 + a + b - c)^2;
Jsdmsolowcd2(4,3) = (a + b) / (1 + a + b - c)^2;
Vsdmsolowcd2 = sdmsolowcd2.minfo(2:4,2:4);
varsdmsolowcd2 = Jsdmsolowcd2*Vsdmsolowcd2*Jsdmsolowcd2';
sigalphasdmsolowcd2 = sqrt(varsdmsolowcd2(1,1));
sigphisdmsolowcd2 = sqrt(varsdmsolowcd2(2,2));
siggammasdmsolowcd2 = sqrt(varsdmsolowcd2(3,3));
sigREsdmsolowcd2 = sqrt(varsdmsolowcd2(4,4));
talphasdmsolowcd2 = abs(alphasdmsolowcd2/sigalphasdmsolowcd2);
tphisdmsolowcd2 = abs(phisdmsolowcd2/sigphisdmsolowcd2);
tgammasdmsolowcd2 = abs(gammasdmsolowcd2/siggammasdmsolowcd2);
tREsdmsolowcd2 = abs(REsdmsolowcd2/sigREsdmsolowcd2);
palphasdmsolowcd2 = 2*(1 - norm_cdf (talphasdmsolowcd2))
pphisdmsolowcd2 = 2*(1 - norm_cdf (tphisdmsolowcd2))
pgammasdmsolowcd2 = 2*(1 - norm_cdf (tgammasdmsolowcd2))
pREsdmsolowcd2 = 2*(1 - norm_cdf (tREsdmsolowcd2))

% Structural parameters (We2)

constant = sdmsolowce2.beta(1);
a = sdmsolowce2.beta(2);
b = sdmsolowce2.beta(3);
c = sdmsolowce2.rho;

alphasdmsolowce2 = b / (b - c)
phisdmsolowce2 = a / (1 + a) - alphasdmsolowce2
gammasdmsolowce2 = - (b - c) / (1 + a)
REsdmsolowce2 = (a + b) / (1 + a + b - c)
Jsdmsolowce2=zeros(4,3);
Jsdmsolowce2(1,1) = 0;
Jsdmsolowce2(1,2) = -c / (b - c)^2;
Jsdmsolowce2(1,3) = b / (b - c)^2;
Jsdmsolowce2(2,1) = 1 / (1 + a)^2;
Jsdmsolowce2(2,2) = c / (b - c)^2;
Jsdmsolowce2(2,3) = - b / (b - c)^2;
Jsdmsolowce2(3,1) = (b - c) / (1 + a)^2;
Jsdmsolowce2(3,2) = - 1 / (1 + a);
Jsdmsolowce2(3,3) = 1 / (1 + a);
Jsdmsolowce2(4,1) = (1 - c) / (1 + a + b - c)^2;
Jsdmsolowce2(4,2) = (1 - c) / (1 + a + b - c)^2;
Jsdmsolowce2(4,3) = (a + b) / (1 + a + b - c)^2;
Vsdmsolowce2 = sdmsolowce2.minfo(2:4,2:4);
varsdmsolowce2 = Jsdmsolowce2*Vsdmsolowce2*Jsdmsolowce2';
sigalphasdmsolowce2 = sqrt(varsdmsolowce2(1,1));
sigphisdmsolowce2 = sqrt(varsdmsolowce2(2,2));
siggammasdmsolowce2 = sqrt(varsdmsolowce2(3,3));
sigREsdmsolowce2 = sqrt(varsdmsolowce2(4,4));
talphasdmsolowce2 = abs(alphasdmsolowce2/sigalphasdmsolowce2);
tphisdmsolowce2 = abs(phisdmsolowce2/sigphisdmsolowce2);
tgammasdmsolowce2 = abs(gammasdmsolowce2/siggammasdmsolowce2);
tREsdmsolowce2 = abs(REsdmsolowce2/sigREsdmsolowce2);
palphasdmsolowce2 = 2*(1 - norm_cdf (talphasdmsolowce2))
pphisdmsolowce2 = 2*(1 - norm_cdf (tphisdmsolowce2))
pgammasdmsolowce2 = 2*(1 - norm_cdf (tgammasdmsolowce2))
pREsdmsolowce2 = 2*(1 - norm_cdf (tREsdmsolowce2))

% LR test for theoretical constraint

logliksdmsolowd2 = sdmsolowd2.lik;
logliksdmsolowcd2 = sdmsolowcd2.lik;
LRsdmsolowd2 = -2*(logliksdmsolowcd2 - logliksdmsolowd2)
psdmsolowd2 = 1-chis_prb(LRsdmsolowd2,2)

logliksdmsolowe2 = sdmsolowe2.lik;
logliksdmsolowce2 = sdmsolowce2.lik;
LRsdmsolowe2 = -2*(logliksdmsolowce2 - logliksdmsolowe2)
psdmsolowe2 = 1-chis_prb(LRsdmsolowe2,2)

% LR test for common factor

logliksemsolowcd2 = semsolowcd2.lik;
logliksdmsolowcd2 = sdmsolowcd2.lik;
LRsolowcd2 = -2*(logliksemsolowcd2 - logliksdmsolowcd2)
psolowcd2 = 1-chis_prb(LRsolowcd2,1)

logliksemsolowce2 = semsolowce2.lik;
logliksdmsolowce2 = sdmsolowce2.lik;
LRsolowce2 = -2*(logliksemsolowce2 - logliksdmsolowce2)
psolowce2 = 1-chis_prb(LRsolowce2,1)

%%%%%%%%%%%%%%%%%%%%%% Convergence Solow SDM unconstrained

xsolowconv = [ones(n,1) lny60 lns lnngd];
vsolowconv = strvcat('gy','constant','lny60','lns','lnngd');

sdmsolowconvd2 = sdm(gy,xsolowconv,Wd2,info);
prt(sdmsolowconvd2,vsolowconv);
sdmsolowconve2 = sdm(gy,xsolowconv,We2,info);
prt(sdmsolowconve2,vsolowconv);

% Structural parameters (Wd2)

constant = sdmsolowconvd2.beta(1);
a = sdmsolowconvd2.beta(2);
sigsdmsolowconvd2 = sdmsolowconvd2.sige;
t = 35;
lambdasdmsolowconvd2 = - (1/t) * log(1+a*t)
Jsdmsolowconvd2 = - 1 / (1+a*t);
sigsdmsolowconvd2 = abs(sigsdmsolowconvd2*Jsdmsolowconvd2);
tlambdasdmsolowconvd2 = abs(lambdasdmsolowconvd2 / sigsdmsolowconvd2);
plambdasdmsolowconvd2 = 2*(1 - norm_cdf (tlambdasdmsolowconvd2))

% Structural parameters (We2)

constant = sdmsolowconve2.beta(1);
a = sdmsolowconve2.beta(2);
sigsdmsolowconve2 = sdmsolowconve2.sige;
t = 35;
lambdasdmsolowconve2 = - (1/t) * log(1+a*t)
Jsdmsolowconve2 = - 1 / (1+a*t);
sigsdmsolowconve2 = abs(sigsdmsolowconve2*Jsdmsolowconve2);
tlambdasdmsolowconve2 = abs(lambdasdmsolowconve2 / sigsdmsolowconve2);
plambdasdmsolowconve2 = 2*(1 - norm_cdf (tlambdasdmsolowconve2))

% LR test for common factor

logliksemsolowconvd2 = semsolowconvd2.lik;
logliksdmsolowconvd2 = sdmsolowconvd2.lik;
LRsolowconvd2 = -2*(logliksemsolowconvd2 - logliksdmsolowconvd2)
psolowconvd2 = 1-chis_prb(LRsolowconvd2,3)

logliksemsolowconve2 = semsolowconve2.lik;
logliksdmsolowconve2 = sdmsolowconve2.lik;
LRsolowconve2 = -2*(logliksemsolowconve2 - logliksdmsolowconve2)
psolowconve2 = 1-chis_prb(LRsolowconve2,3)

