% Developed by J.Paul Elhorst 2011
% University of Groningen
% Department of Economics, Econometrics and Finance
% P.O. Box 800
% 9700AV Groningen
% the Netherlands
% j.p.elhorst@rug.nl
%
% REFERENCE: 
% Elhorst J.P., Lacombe D.J., Piras G. (2012) 
% On Model Specification and Parameter Space Definitions in Higher Order Spatial Econometrics Models. 
% Regional Science and Urban Economics 42: 211-220.

% Determination of W matrices
% All based on three-dimensional Bucky Ball with N=60
% Due to symmetry, question whether to row-normalize or to normalize by
% largest eigenvalue does not affect the results
%
N=60;
% Binary Contiguity matrices
WBC1=normw(bucky); % Binary contiguity matrix of first-order neighbors 
WBC2=normw(bucky*bucky-3*eye(N)); % Binary contiguity matrix of second-order neighbors only
% Nearest Neighbor matrices
WNN3=WBC1;         % Also 3 nearest neighbor matrix
WNN9=normw(bucky*bucky-3*eye(N)+bucky); % 9 nearest neighbor matrix
% Group interaction matrices
WG1=full(bucky);
WG2=full(bucky);
WG1(41:60,41:60)=zeros(20,20);
WG1(1:40,41:60)=zeros(40,20);
WG1(41:60,1:40)=zeros(20,40);
WG2(1:40,1:40)=zeros(40,40);
WG2(1:40,41:60)=zeros(40,20);
WG2(41:60,1:40)=zeros(20,40);
WG1=normw(WG1);
WG2=normw(WG2);
% Construction of inverse distance matrix
W=full(bucky);
Wp=W;
p=min(Wp(:,1));
count=1;
while p==0
    Wp=W*Wp;
    count=count+1;
    p=min(Wp(:,1));
end
count
WID=zeros(N,N);
WID2=zeros(N,N);
Wsign=zeros(N,N);
WBC3=zeros(N,N);
for c=1:count
    if (c==1) Wlast=zeros(N,N);Wnext=W; 
    else Wlast=Wnext; Wnext=W*Wnext;
    end
for i=1:N
    for j=1:N
        if (Wlast(i,j)~=Wnext(i,j) && Wsign(i,j)==0 && i~=j)
            WID(i,j)=1/c;
            Wsign(i,j)=1;
        end
    end
end
end

% Determination WBC3 matrix
for i=1:N
    for j=1:N
        if (WID(i,j)==1/3) WBC3(i,j)=1; end
    end
end
WBC3=normw(WBC3);

% Determination WID2 matrix
for i=1:N
    for j=1:N
        if (WID(i,j)>1/2.5) WID2(i,j)=WID(i,j); end
    end
end

WID=normw(WID);
WID2=normw(WID2);
%
% Smallest eigenvalues and their reciprocals
%
[min(eig(WBC1)) 1/min(eig(WBC1))];
[min(eig(WBC2)) 1/min(eig(WBC2))]
[min(eig(WBC3)) 1/min(eig(WBC3))]
[min(eig(WNN3)) 1/min(eig(WNN3))]
[min(eig(WNN9)) 1/min(eig(WNN9))]
[min(eig(WG1)) 1/min(eig(WG1))]
[min(eig(WG2)) 1/min(eig(WG2))]
[min(eig(WID)) 1/min(eig(WID))]
[min(eig(WID2)) 1/min(eig(WID2))]
%
% Naive figure, using deltaone
% Time-series form, using deltatime
%
deltaone=zeros(360,2);
deltatime=zeros(360,2);
for a=1:90
    deltaone(a,1)=1-a/90;
    deltaone(a,2)=-a/90;
    if (a<30) 
        deltatime(a,1)=1+a/30;
        deltatime(a,2)=-a/30;  
    else
        deltatime(a,1)=2-(a-30)/30;
        deltatime(a,2)=-1;  
    end
end
for a=91:180
    deltaone(a,1)=-(a-90)/90;
    deltaone(a,2)=-1+(a-90)/90;
    if (a<150) 
        deltatime(a,1)=-(a-90)/30;
        deltatime(a,2)=-1;  
    else
        deltatime(a,1)=-2+(a-150)/30;
        deltatime(a,2)=-1+(a-150)/30;  
    end
end    
for a=181:270
    deltaone(a,1)=-1+(a-180)/90;
    deltaone(a,2)=(a-180)/90;
    deltatime(a,1)=-1+(a-180)/90;
    deltatime(a,2)=(a-180)/90;  
end    
for a=270:360
    deltaone(a,1)=0+(a-270)/90;
    deltaone(a,2)=1-deltaone(a,1);
    deltatime(a,1)=(a-270)/90;
    deltatime(a,2)=1-(a-270)/90;  
end
%
% Combinations of TWO different W matrices
%
% Two binary contiguity matrices, no overlap
%
% Define W1 and W2
W1=WBC1;
W2=WBC2;
sum(sum(W1.*W2))
delta12=zeros(360,2);
% second quadrant
for a=1:89
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
a=90; % Note tand(90) is not defined
    eigmin=min(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmin;
% Third quadrant
for a=91:180
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% Fourth quadrant
for a=181:269
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% First quadrant
a=270; % Note tand(270) is not defined
    eigmax=max(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmax;
for a=271:360
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
figure;
plot(delta12(:,1),delta12(:,2),'b*',deltatime(:,1),deltatime(:,2),'g+',deltaone(:,1),deltaone(:,2),'r.');
grid on;
%
% Two group interaction matrices, no overlap
% WG1 and WG2 also satisfy Craig-Sakamoto lemma,
% which gives a rectangle
%
W1=WG1;
W2=WG2;
sum(sum(W1.*W2))
delta12=zeros(360,2);
% second quadrant
for a=1:89
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
a=90;
    eigmin=min(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmin;
% Third quadrant
for a=91:180
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% Fourth quadrant
for a=181:269
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% First quadrant
a=270;
    eigmax=max(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmax;
for a=271:360
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
figure;
plot(delta12(:,1),delta12(:,2),'b*',deltatime(:,1),deltatime(:,2),'g+',deltaone(:,1),deltaone(:,2),'r.');
grid on;
%
% Matrices that overlap
%
% 3 and 9 nearest neighbor matrices
%
% 3 nearest neighbor matrix = WBC1
% WBC2=WNN9-WNN3
%
W1=WNN3;
W2=WNN9;
sum(sum(W1.*W2))
delta12=zeros(360,2);
% second quadrant
for a=1:89
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
a=90;
    eigmin=min(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmin;
% Third quadrant
for a=91:180
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% Fourth quadrant
for a=181:269
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% First quadrant
a=270;
    eigmax=max(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmax;
for a=271:360
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
figure;
plot(delta12(:,1),delta12(:,2),'b*',deltatime(:,1),deltatime(:,2),'g+',deltaone(:,1),deltaone(:,2),'r.');
grid on;
%
% WNN9 and WID2
%
W1=WNN9;
W2=WID2;
sum(sum(W1.*W2))
delta12=zeros(360,2);
% second quadrant
for a=1:89
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
a=90;
    eigmin=min(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmin;
% Third quadrant
for a=91:180
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% Fourth quadrant
for a=181:269
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% First quadrant
a=270;
    eigmax=max(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmax;
for a=271:360
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
figure;
plot(delta12(:,1),delta12(:,2),'b*',deltatime(:,1),deltatime(:,2),'g+',deltaone(:,1),deltaone(:,2),'r.');
grid on;
%
% WID2 and WID
%
W1=WID2;
W2=WID;
sum(sum(W1.*W2))
delta12=zeros(360,2);
% second quadrant
for a=1:89
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
a=90;
    eigmin=min(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmin;
% Third quadrant
for a=91:180
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% Fourth quadrant
for a=181:269
    eigmin=min(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmin;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
% First quadrant
a=270;
    eigmax=max(real(eig(W2)));
    delta12(a,1)=0;
    delta12(a,2)=1/eigmax;
for a=271:360
    eigmax=max(real(eig(W1+tand(-a)*W2)));
    delta12(a,1)=1/eigmax;
    delta12(a,2)=tand(-a)*delta12(a,1);
end
figure;
plot(delta12(:,1),delta12(:,2),'b*',deltatime(:,1),deltatime(:,2),'g+',deltaone(:,1),deltaone(:,2),'r.');
grid on;
%
%
% Direct/Indirect effects estimates
%
%
rho1=-0.2;rho2=0.9;
S=inv(eye(N)-rho1*WNN3-rho2*WNN9);
EAVD=sum(diag(S))/N; % average direct effect
EAVI=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC=sum(sum(S,1)'-diag(S))/N; % average indirect effect
[EAVD EAVI EAVC]
rho1=0.1;rho2=0.6;
S=inv(eye(N)-rho1*WBC1-rho2*WBC2);
EAVD=sum(diag(S))/N; % average direct effect
EAVI=sum(sum(S,2)-diag(S))/N; % average indirect effect
EAVC=sum(sum(S,1)'-diag(S))/N; % average indirect effect
[EAVD EAVI EAVC]
% Magnitude order-effects
IS=rho1*WBC1+rho2*WBC2;
sum(sum(IS,2)-diag(IS))/N
IS2=IS*IS;
sum(sum(IS2,2)-diag(IS2))/N
IS3=IS*IS2;
sum(sum(IS3,2)-diag(IS3))/N
%
% Second-order polynomial
%
SP=inv(eye(N)-rho2*WBC2)*inv(eye(N)-rho1*WBC1);
EAVD=sum(diag(SP))/N; % average direct effect
EAVI=sum(sum(SP,2)-diag(SP))/N; % average indirect effect
EAVC=sum(sum(SP,1)'-diag(SP))/N; % average indirect effect
[EAVD EAVI EAVC]
%
SP10=inv(eye(N)-rho1*WBC2)-eye(N);
sum(sum(SP10,2)-diag(SP10))/N
SP11=rho1*WBC1;
sum(sum(SP11,2)-diag(SP11))/N
SP12=SP11*(rho1*WBC1);
sum(sum(SP12,2)-diag(SP12))/N
SP13=SP12*(rho1*WBC1);
sum(sum(SP13,2)-diag(SP13))/N
%
SP20=inv(eye(N)-rho2*WBC2)-eye(N);
sum(sum(SP20,2)-diag(SP20))/N
SP21=rho2*WBC2;
sum(sum(SP21,2)-diag(SP21))/N
SP22=SP21*(rho2*WBC2);
sum(sum(SP22,2)-diag(SP22))/N
SP23=SP22*(rho2*WBC2);
sum(sum(SP23,2)-diag(SP23))/N
%
% Combinations of THREE different W matrices
%
% Three binary contiguity matrices, no overlap
%
% Define W1, W2 and W3, delta3 is held constant
%
delta3four=[-0.6;-0.2;0.2;0.6]; % four values of delta 3 are considered
W1=WBC1;
W2=WBC2;
W3=WBC3;
delta123=zeros(4,360,3);
deltaone=zeros(4,360,3); % Naive figure, using deltaone
for k=1:4
%   a2=alpha(k);
    delta3=delta3four(k);
% second quadrant
for a1=1:89
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmax=max(real(eig(W1+tand(-a1)*W2+tand(a2)/cosd(-a1)*W3)));
        delta123(k,a1,1)=1/eigmax;
        delta123(k,a1,2)=tand(-a1)*delta123(k,a1,1);
        delta123(k,a1,3)=tand(a2)/cosd(-a1)*delta123(k,a1,1);
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end
    deltaone(k,a1,1)=(1-a1/90)*(1-abs(delta3));
    deltaone(k,a1,2)=(-a1/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
end
a1=90; % Note tand(90) is not defined
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmin=min(real(eig(W2-tand(a2)*W3)));
        delta123(k,a1,1)=0;
        delta123(k,a1,2)=1/eigmin;
        delta123(k,a1,3)=-tand(a2)*delta123(k,a1,2);    
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end
    deltaone(k,a1,1)=(1-a1/90)*(1-abs(delta3));
    deltaone(k,a1,2)=(-a1/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
% Third quadrant
for a1=91:180
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmin=min(real(eig(W1+tand(-a1)*W2+tand(a2)/cosd(-a1)*W3)));
        delta123(k,a1,1)=1/eigmin;
        delta123(k,a1,2)=tand(-a1)*delta123(k,a1,1);
        delta123(k,a1,3)=tand(a2)/cosd(-a1)*delta123(k,a1,1);
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end
    deltaone(k,a1,1)=(-(a1-90)/90)*(1-abs(delta3));
    deltaone(k,a1,2)=(-1+(a1-90)/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
end
% Fourth quadrant
for a1=181:269
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmin=min(real(eig(W1+tand(-a1)*W2+tand(a2)/cosd(-a1)*W3)));
        delta123(k,a1,1)=1/eigmin;
        delta123(k,a1,2)=tand(-a1)*delta123(k,a1,1);
        delta123(k,a1,3)=tand(a2)/cosd(-a1)*delta123(k,a1,1);
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end   
    deltaone(k,a1,1)=(-1+(a1-180)/90)*(1-abs(delta3));
    deltaone(k,a1,2)=((a1-180)/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
end
% First quadrant
a1=270; % Note tand(270) is not defined
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmax=max(real(eig(W2+tand(a2)*W3)));
        delta123(k,a1,1)=0;
        delta123(k,a1,2)=1/eigmax;
        delta123(k,a1,3)=tand(a2)*delta123(k,a1,2);    
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end
    deltaone(k,a1,1)=(-1+(a1-180)/90)*(1-abs(delta3));
    deltaone(k,a1,2)=((a1-180)/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
for a1=271:360
    if (delta3<0) 
        x0=-0.001;
        x1=-89.999;
    else
        x0=0.001;
        x1=89.999;
    end    
    precision=0.0005;
    fprecision=1000;
    count=0;
    while fprecision>precision
        count=count+1;
        a2=0.5*x0+0.5*x1;
        eigmax=max(real(eig(W1+tand(-a1)*W2+tand(a2)/cosd(-a1)*W3)));
        delta123(k,a1,1)=1/eigmax;
        delta123(k,a1,2)=tand(-a1)*delta123(k,a1,1);
        delta123(k,a1,3)=tand(a2)/cosd(-a1)*delta123(k,a1,1);
        f=delta123(k,a1,3)-delta3;
        if (f>0 && delta3<0) x0=a2; 
        elseif (f<0 && delta3<0) x1=a2; 
        elseif (f>0 && delta3>0) x1=a2; 
        else x0=a2; 
        end
        fprecision=abs(f);
    end
    deltaone(k,a1,1)=((a1-270)/90)*(1-abs(delta3));
    deltaone(k,a1,2)=(1-(a1-270)/90)*(1-abs(delta3));
    deltaone(k,a1,3)=delta3;
end
delta(:,:)=delta123(k,:,:);
deltao(:,:)=deltaone(k,:,:);
figure;
plot(delta(:,1),delta(:,2),'b*',deltao(:,1),deltao(:,2),'r.');
grid on;
end