% This Program Checks the Proofs in the Paper Sticky Prices in a Dynamic Network Economy

% This program checks the proofs in the paper "Sticky Prices in a Dynamic Network Economy:

% A Family of Counterexamples."

%

% Propositions 1,2 and 4 are striaghtforward and can be checked for correctness

% in the usual fashion. They are not handled by this program. Most other propositions

% in the paper assert that a particuar equality is true. For instance,

%

% Claim: LHS(x1,x2) = RHS(x1,x2)

%

% where x1 and x2 are parameters. Clearly the details of LHS and RHS depend on

% the context. However, in most cases, LHS and RHS can be computed in different

% ways (otherwise there would be nothing to prove). Therefore, to check the claim,

% it suffices to randomly choose values for x1 and x2, and then compute LHS(x1,x2)

% one way, and RHS(x1,x2) in another way. If the two computations agree up to 14 decimal

% places, for any choices of x1 and x2, then it is almost impossible for the claim to be

% false. (Of course, I suppose that this could also mean that the claim is

% trivially true, but that is not really the problem I am worried about.)

%

% **** THE PROGRAM CHECKS THE SPECIAL CASE OF F={z1,z2} ... THAT IS, THERE ARE EXACTLY TWO

% **** STICKY PRICES.

%

%

% The program makes several calls to a simple function called "ceck_matrix( ),"

% which takes as input an m x n matrix. This function needs to be stored

% in a place where it can be accessed by this program. The function is

% printed just below (for cutting and pasting), and it can also be found at

%

%

%

%

% HERE IS checkmatrix(A)

%

%function E=checkmatrix(A)

%%Checks whether any of the entries of A exceeds 10^-12

%

%[m,n]=size(A);

%

%cc=.5*ones(m,n);

%dd=zeros(m,n);

%for i=1:m

% for j=1:n

% if abs(A(i,j))>10^-12

% dd(i,j)=1;

% else

% dd(i,j)=0;

% end

% end

%end

% tt=(dd>cc);

% if any(tt)

% E=1;

% else

% E=0;

% end

%

% ***************************************************************

% ***************************************************************

%

% The program can be run in a normal Matlab session. To run the program, first assign the

% following variables:

%

% p - a prime number equal to 1 mod 4 (this is the number of vertices in graph)

% a - a real number greater than 1 (corresponds to $\alpha$ in the paper).

% x - any number between -1 and 1 satisfying |x|/(2*(1-x))<1 (so-2/3<x<2/3 will do)

% t - a positive integer (a time period)

% t1 - an odd positive integer (an odd-numbered time period)

% t2 - an even positive integer (an even-numbered time period)

% z1 and z2 - elements of the set {1,2,...,p-1}. Moreover z1 and z2 must be

% such that z1-z2 is NOT be a non-zero square mod p (that is, z1-z2 cannot

% be a quadratic residue mod p), and both z1 and z2 must either belong to N(0),

% or they must not belong to N[0]. (That is, EITHER z1 and z2 are both

% quadratic residues, OR neither is a quadratic residues.)

%

%

% Now run: DNE_sticky_check.

%

% Here is an example:

%

% > p=29; a=2.3; x=.4572197;

% > t=24; t1=11;t2=8;

% > z1=2;z2=10;

% > p_fix_check2

%

% A succesful run looks like this

% Step 1 okay.

% Step 2.1 okay.

% Step 2.2 okay.

% Constants in Step 2.3 okay.

% Step 2.3 okay.

% Step 2.4 okay.

% Step 2.5 okay.

% Step 3 okay.

% mu okay.

% Step 4 okay.

% eta okay.

% Step 5 okay.

% Step 6 okay.

% Step 7 okay for odd t.

% Step 7 okay for even t.

% Step 8 okay.

% Here is what the steps do:

% Step 1 checks validity of eq'n (21)

% Step 2, Par ts 2.1, 2.2, 2.3, 2.4, and 2.5 check Parts 1 - 5 of Lemma 5

% Step 3 checks Lemma 6

% Step 4 checks the values of mu1 and mu2 (mu is calculated in 3 ways: (i). Part 3

% of Lemma 5, (ii). eq'n (35), and (iii). according to polynomials in Section 7.4.4.)

% Step 5 checks the values of eta1 and eta2 (eta is calculated in 3 ways: (i). Part 3

% of Lemma 5, (ii). eq'n (34), and (iii). according to polynomials in Section 7.4.4.)

% Step 6 checks Lemma 8

% Step 7 checks Proposition 3. In Case 1 t is odd, and in Case 2 t is even.

% Step 8 is redundant. It checks a computation that has been removed from the paper.

% CHECKMATRIX

% The program makes frequent calls to the function "checkmatrix," which

% follows. You will need to remove the function from this file,

% and store it where it can be invoked.

% SPATIAL_STICKY_CHECK BEGINS:

% builds the adjacency matrix for

% a Paley graph on p vertices

q1=mod(p,4);

q2=isprime(p);

if q1~=1 | q2~=1

disp('EITHER p IS NOT A PRIME, OR p is NOT equal TO 1 MOD 4')

end

for i=1:p-1

h0(i)=mod(i^2,p);

end

h1=sort(h0);

for i=1:(p-1)/2

s(i)=h1(2*i);

end

c1=zeros(p,1);

r=1;

for i=1:p-1

if i==s(r)

c1(i+1,1)=1;

r=r+1;

else

c1(i+1,1)=0;

end

end

c_next=c1;

P_adj=[c_next];

for i=2:p

c_hold=c_next;

c_next(1,1)=c_hold(p,1);

for i=1:p-1

c_next(i+1,1)=c_hold(i,1);

end

P_adj=[P_adj c_next];

end

%Checks whether z1 and z2 are neighbours

for i=1:p

if abs(z1-z2)==mod(i^2,p)

tst(i)=1;

else

tst(i)=0;

end

end

if sum(tst)>.5

disp('z1 and z2 are neighbours: program will likely fail')

end

if P_adj(z1+1,1)-P_adj(z2+1,1)~=0

disp('Either z1 is a neighbour of 0 and z2 is not, or vice versa; program will likely fail')

end

%Basic parameters

A=1/((p-1)/2)*P_adj;

%a=2;

s=.2*p;

ms=(p+s)/p;

%Eigenvalues and eigenvectors

evl0=1;

evl1=(sqrt(p)-1)/(p-1);

evl2=(-sqrt(p)-1)/(p-1);

w0=(p-1)/2;

w1=1/2*(sqrt(p)-1);

w2=1/2*(-sqrt(p)-1);

%The matrrices (these have the same names as in the paper).

P1=zeros(p,p);

P2=zeros(p,p);

P1(z1+1,z1+1)=1;

P2(z2+1,z2+1)=1;

P=P1+P2;

B=A-A*P+1/a*A*P*A;

C=eye(p)+1/a*A*P;

U1=-A*P1*A+1/a*A*P1*A^2;

U2=-A*P2*A+1/a*A*P2*A^2;

U=-A*P*A+1/a*A*P*A^2;

E=eye(p)-P+1/a*P*A;

D=E*A;

K1=inv(eye(p)-(evl1/a)^2*A^2);

L1=inv(eye(p)-(evl1/a)^2*B*A);

K2=inv(eye(p)-(evl2/a)^2*A^2);

L2=inv(eye(p)-(evl2/a)^2*B*A);

M1=inv(eye(p)-(evl1/a)^2*D*A);

M2=inv(eye(p)-(evl2/a)^2*D*A);

Kx=inv(eye(p)-x*A^2); %K and L from the proof of Lemma 5

Lx=inv(eye(p)-x*B*A);

Nz1=P_adj(:,z1+1); %neighbours of z1 \in {0,...,p-1}

Nz2=P_adj(:,z2+1); %neighbours of z2 \in {0,...,p-1}

Nz=Nz1+Nz2;

F1=zeros(p,1); %the vector 1_z

F2=zeros(p,1); %the vector 1_z

F1(z1+1,1)=1;

F2(z2+1,1)=1;

F=F1+F2;

e1=zeros(p,1);

e1(1)=1;

N0=P_adj(:,1); %the vecto 1_N(0)

D0=ones(p,1)-N0-e1; %the vector 1_N[0]

v1=w0*e1+w1*N0+w2*D0; %v_1 and v_2 in eq'n (21)

v2=w0*e1+w2*N0+w1*D0;

m1_ppr=ms*(ones(p,1)+s/(p+s)*(v1+v2)); %m_1 in (22)

m1=ones(p,1);

m1(1)=1+s; %m_1 from (9)

%**********************************************************************

%Checking

%**********************************************************************

eig_mlb=eig(A); %eigenvalues of A computed by MATLAB compare with evl0, evl1 and evl2

% on last line of section 7.1.2

%**********************************************************************

%Step 1: Check validity of eq'n (21)

%**********************************************************************

lhs1=A^t*m1; %Computes m_{t+1} as in (10).

rhs1=ms*(ones(p,1)+s/(p+s)*(evl1^t*v1+evl2^t*v2)); %Computes m_{t+1} as in (22)

diff1=lhs1-rhs1;

check_error=checkmatrix(diff1);

if check_error>.5

disp('Error in Step 1.')

else

disp('Step 1 okay.')

end

%**********************************************************************

%Step 2: Check Lemma 5

%**********************************************************************

%*********************

%Part 1

lhs2a=Kx-Lx; %Computes LHS of Lemma 5, Part 1.

rhs2a=-x*Kx*U*Lx; %Computes RHS of Lemma 5, Part 1.

diff2a=lhs2a-rhs2a;

check_error=checkmatrix(diff2a);

if check_error>.5

disp('Error in Step 2.1')

else

disp('Step 2.1 okay.')

end

%*********************

%Part 2

rhs2b=rank(rhs2a); %Computes rank of K-L.

if rhs2b==2

disp('Step 2.2 okay.')

else

disp('Error in Step 2.2')

end

%*********************

%Part 3

%Computes rhs of eq'n in Part 3: All combinations.

lhsc1=U1*Kx*Nz1;

lhsc2=U1*Kx*Nz2;

lhsc3=U2*Kx*Nz2;

lhsc4=U2*Kx*Nz1;

lhsc5=U*Kx*Nz;

const1=sum(U1*Kx*Nz1)/sum(Nz1);

const2=sum(U1*Kx*Nz2)/sum(Nz1);

const3=sum(U2*Kx*Nz2)/sum(Nz2);

const4=sum(U2*Kx*Nz1)/sum(Nz2);

const5=sum(U*Kx*Nz)/sum(Nz);

const6=const1+const2;

%Computes rhs's

rhsc1=const1*Nz1;

rhsc2=const2*Nz1;

rhsc3=const3*Nz2;

rhsc4=const4*Nz2;

rhsc5=const5*Nz;

diffc1=lhsc1-rhsc1;

diffc2=lhsc2-rhsc2;

diffc3=lhsc3-rhsc3;

diffc4=lhsc4-rhsc4;

diffc5=lhsc5-rhsc5;

de1=(const1-const3>10^-12);

de2=(const2-const4>10^-12);

de3=(const5-const6>10^-12);

if de1+de2+de3>.5

disp('Problem with constants in Step 2.3')

else

disp('Constants in Step 2.3 okay.')

end

ce1=checkmatrix(diffc1);

ce2=checkmatrix(diffc2);

ce3=checkmatrix(diffc3);

ce4=checkmatrix(diffc4);

ce5=checkmatrix(diffc5);

if ce1+ce2+ce3+ce4+ce5>.5

disp('Error in Step 2.3')

else

disp('Step 2.3 okay.')

end

%*********************

%Part 4

lhsd=Lx*Nz; %Computes lhs of eq'n in Part 4 for case |F|=2.

const=1/(1-x*(const1+const2)); %computes const on rhs; const1 and const2 (which correspond to mu and eta) are retrieved from previous step.

rhsd=const*Kx*Nz;

diffd=lhsd-rhsd;

check_error=checkmatrix(diffd);

if check_error>.5

disp('Error in Step 2.4')

else

disp('Step 2.4 okay.')

end

%*********************

%Part 5

lhse1=(Kx-Lx)*v1; %Computes lhs of eq'n in Part 5, for case |F|=2.

const=sum((Kx-Lx)*v1)/sum(Kx*Nz); %Computes rhs.

rhse1=const*Kx*Nz;

lhse2=(Kx-Lx)*v2;

const=sum((Kx-Lx)*v2)/sum(Kx*Nz);

rhse2=const*Kx*Nz;

diffe1=lhse1-rhse1;

diffe2=lhse2-rhse2;

check_error1=checkmatrix(diffe1);

check_error2=checkmatrix(diffe2);

if check_error1+check_error2>.5

disp('Error in Step 2.5')

else

disp('Step 2.5 okay.')

end

%**********************************************************************

%Step 3: Check Lemma 6

%**********************************************************************

if P_adj(z1+1,1)==1 & P_adj(z2+1,1)==1

lhs3a=A*P*v1;

rhs3a=evl1*Nz;

lhs3b=A*P*v2;

rhs3b=evl2*Nz;

else

lhs3a=A*P*v1;

rhs3a=evl2*Nz;

lhs3b=A*P*v2;

rhs3b=evl1*Nz;

end

diff3a=lhs3a-rhs3a;

diff3b=lhs3b-rhs3b;

check_error1=checkmatrix(diff3a);

check_error2=checkmatrix(diff3a);

if check_error1+check_error2>.5

disp('Error in Step 3')

else

disp('Step 3 okay.')

end

%**********************************************************************

%Step 4: Check mu

%**********************************************************************

%the entries of K(x) from Section 7.4.1

x1=evl1^2/a^2;

kz1=1/p*(1/(1-x1)+(p-1)/2*1/(1-x1*evl1^2)+(p-1)/2*1/(1-x1*evl2^2));

ka1=1/p*(1/(1-x1)+(sqrt(p)-1)/2*1/(1-x1*evl1^2)+(-sqrt(p)-1)/2*1/(1-x1*evl2^2));

kn1=1/p*(1/(1-x1)+(-sqrt(p)-1)/2*1/(1-x1*evl1^2)+(sqrt(p)-1)/2*1/(1-x1*evl2^2));

x2=evl2^2/a^2;

kz2=1/p*(1/(1-x2)+(p-1)/2*1/(1-x2*evl1^2)+(p-1)/2*1/(1-x2*evl2^2));

ka2=1/p*(1/(1-x2)+(sqrt(p)-1)/2*1/(1-x2*evl1^2)+(-sqrt(p)-1)/2*1/(1-x2*evl2^2));

kn2=1/p*(1/(1-x2)+(-sqrt(p)-1)/2*1/(1-x2*evl1^2)+(sqrt(p)-1)/2*1/(1-x2*evl2^2));

col1_1=kz1*e1+ka1*N0+kn1*D0; %Column 1 of K1(x)

col1_2=kz2*e1+ka2*N0+kn2*D0; %Column 1 of K2(x)

%uz, ua, and un from Section 7.4.1

uz=4/a*1/(p-1)^2;

ua=2/a*(p-5)/(p-1)^3-4/(p-1)^2;

un=2/a*1/(p-1)^2;

%the numbers A_/rho, B_/rho, and C_/rho, from Section 7.4.2

A1=(p-1)/2*ka1;

B1=kz1+(p-5)/4*ka1+(p-1)/4*kn1;

C1=(p-1)/4*ka1+(p-1)/4*kn1;

A2=(p-1)/2*ka2;

B2=kz2+(p-5)/4*ka2+(p-1)/4*kn2;

C2=(p-1)/4*ka2+(p-1)/4*kn2;

%Compute mu from 7.4.3 (eq'n (35))

mu1=A1*uz+(p-1)/2*B1*ua+(p-1)/2*C1*un;

mu2=A2*uz+(p-1)/2*B2*ua+(p-1)/2*C2*un;

%Computes mu as in Section 7.4.4 (e.g mu1 is from eq'n 36).

n_mu1=2*a^4*(sqrt(p)+1)^4*(p-1)-a^3*(sqrt(p)+1)^4*(p-5)-a^2*(sqrt(p)+1)^2*(p+3)-a*(sqrt(p)+1)^2+1;

d1=a^4*(sqrt(p)+1)^4*(p-1)^2-2*a^2*(sqrt(p)+1)^2*(p+1)+1;

n_mu2=2*a^4*(sqrt(p)-1)^4*(p-1)-a^3*(sqrt(p)-1)^4*(p-5)-a^2*(sqrt(p)-1)^2*(p+3)-a*(sqrt(p)-1)^2+1;

d2=a^4*(sqrt(p)-1)^4*(p-1)^2-2*a^2*(sqrt(p)-1)^2*(p+1)+1;

mu1_polyn=-a^2*(sqrt(p)+1)^2/(a^2*(sqrt(p)+1)^2-1)*(n_mu1/d1);

mu2_polyn=-a^2*(sqrt(p)-1)^2/(a^2*(sqrt(p)-1)^2-1)*(n_mu2/d2);

diff_mu1=mu1-mu1_polyn;

diff_mu2=mu2-mu2_polyn;

%Checking consistency among different computations of mu

cmu1=checkmatrix(diff_mu1);

cmu2=checkmatrix(diff_mu2);

if cmu1+cmu2>.5

disp('Error in computing an mu')

else

disp('mu okay.')

end

%Computes mu according to eq'n in Part 3 of Lemma 6. Checks all possibilities.

lhs_mu1a=U1*K1*Nz1;

rhs_mu1a=mu1*Nz1;

lhs_mu1b=U2*K1*Nz2;

rhs_mu1b=mu1*Nz2;

lhs_mu2a=U1*K2*Nz1;

rhs_mu2a=mu2*Nz1;

lhs_mu2b=U2*K2*Nz2;

rhs_mu2b=mu2*Nz2;

diff_mu1a=lhs_mu1a-rhs_mu1a;

diff_mu1b=lhs_mu1b-rhs_mu1b;

diff_mu2a=lhs_mu2a-rhs_mu2a;

diff_mu2b=lhs_mu2b-rhs_mu2b;

ce1=checkmatrix(diff_mu1a);

ce2=checkmatrix(diff_mu1b);

ce3=checkmatrix(diff_mu2a);

ce4=checkmatrix(diff_mu2b);

if ce1+ce2+ce3+ce4>.5

disp('Error in Step 4')

else

disp('Step 4 okay.')

end

%**********************************************************************

%Step 5: Check eta

%**********************************************************************

%Compute eta from 7.4.3, eq'n (36)

eta1=A1*un+C1*uz+(p-1)/4*B1*ua+(p-1)/4*C1*ua+(p-1)/4*B1*un+(p-5)/4*C1*un;

eta2=A2*un+C2*uz+(p-1)/4*B2*ua+(p-1)/4*C2*ua+(p-1)/4*B2*un+(p-5)/4*C2*un;

%Computes eta as in Section 7.4.4.

n_eta1=a^3*(sqrt(p)+1)^2*(p-1)-a^2*(sqrt(p)+1)^2*(p-3)-3*a+1;

n_eta2=a^3*(sqrt(p)-1)^2*(p-1)-a^2*(sqrt(p)-1)^2*(p-3)-3*a+1;

eta1_polyn=-a^3*(sqrt(p)+1)^4/(a^2*(sqrt(p)+1)^2-1)*n_eta1/d1;

eta2_polyn=-a^3*(sqrt(p)-1)^4/(a^2*(sqrt(p)-1)^2-1)*n_eta2/d2;

diff_eta1=eta1-eta1_polyn;

diff_eta2=eta2-eta2_polyn;

%Checking consistency among different computations of eta

ceta1=checkmatrix(diff_eta1);

ceta2=checkmatrix(diff_eta2);

if ceta1+ceta2>.5

disp('Error in computing an eta')

else

disp('eta okay.')

end

%Computes eta according to eq'n in Part 3 of Lemma 6. Checks all possibilities.

lhs_eta1a=U1*K1*Nz2;

rhs_eta1a=eta1*Nz1;

lhs_eta1b=U2*K1*Nz1;

rhs_eta1b=eta1*Nz2;

lhs_eta2a=U1*K2*Nz2;

rhs_eta2a=eta2*Nz1;

lhs_eta2b=U2*K2*Nz1;

rhs_eta2b=eta2*Nz2;

diff_eta1a=lhs_eta1a-rhs_eta1a;

diff_eta1b=lhs_eta1b-rhs_eta1b;

diff_eta2a=lhs_eta2a-rhs_eta2a;

diff_eta2b=lhs_eta2b-rhs_eta2b;

ce1=checkmatrix(diff_eta1a);

ce2=checkmatrix(diff_eta1b);

ce3=checkmatrix(diff_eta2a);

ce4=checkmatrix(diff_eta2b);

if ce1+ce2+ce3+ce4>.5

disp('Error in Step 5')

else

disp('Step 5 okay.')

end

%**********************************************************************

%Step 6: Check Lemma 8

%**********************************************************************

%Computes nu1 and nu2 according to formulas in statement of Lemma 9

nu1=-a*(sqrt(p)+1)*(a*sqrt(p)+a-1)/(a^2*(sqrt(p)+1)^4-1);

nu2=-a*(sqrt(p)-1)*(a*sqrt(p)-a+1)/(a^2*(sqrt(p)-1)^4-1);

nu1_non=-(sqrt(p)+1)/(sqrt(p)-1)*nu1;

nu2_non=-(sqrt(p)-1)/(sqrt(p)+1)*nu2;

%Verifyinbg that these nu1 and nu2 satusfy Part 5 of Lemma 6.

lhs_nu1=(K1-L1)*v1;

lhs_nu2=(K2-L2)*v2;

x1=evl1^2/a^2;

x2=evl2^2/a^2;

if P_adj(z1+1,1)==1 & P_adj(z2+1,1)==1

rhs_nu1=-x1*nu1/(1-x1*(mu1+eta1))*K1*Nz;

rhs_nu2=-x2*nu2/(1-x2*(mu2+eta2))*K2*Nz;

else

rhs_nu1=-x1*nu1_non/(1-x1*(mu1+eta1))*K1*Nz;

rhs_nu2=-x2*nu2_non/(1-x2*(mu2+eta2))*K2*Nz;

end

diff_nu1=lhs_nu1-rhs_nu1;

diff_nu2=lhs_nu2-rhs_nu2;

check_error1=checkmatrix(diff_nu1);

check_error2=checkmatrix(diff_nu2);

if check_error1+check_error2>.5

disp('Error in Step 6.')

else

disp('Step 6 okay.')

end

%**********************************************************************

%Step 7: Check Proposition 3

%**********************************************************************

%computes psi1 and psi2 for the case m=2

psi1=mu1+eta1;

psi2=mu2+eta2;

x1=evl1^2/a^2;

x2=evl2^2/a^2;

aa1=1/(1-x1*psi1); %introduces abreviations aa1 and aa2.

aa2=1/(1-x2*psi2);

%Case 1: t odd

J1=inv(eye(p)-evl1/a*A);

J2=inv(eye(p)-evl2/a*A);

a1_odd=evl1^t1*J1*v1+evl2^t1*J2*v2; %computes D_sigma(log p_t) from (13), t odd (drops constant term)

a2_odd=evl1^t1*L1*(C+evl1/a*B)*v1+evl2^t1*L2*(C+evl2/a*B)*v2; %computes D_sigma(log p_t^F) from (11)

lhs_odd=a1_odd-a2_odd; %computes D_sigma(log p_t - log p_t^F), t odd

%computes D_sigma(log p_t - log p_t^F), t odd, from cases (i) and (ii) in proof of Prop.3

if P_adj(z1+1,1)==1 & P_adj(z2+1,1)==1

term1=1/(sqrt(p)+1)^t1*(aa1*(-sqrt(p))/(a*(sqrt(p)+1)^2-1))*K1*Nz;

term2=1/(sqrt(p)-1)^t1*(aa2*sqrt(p)/(a*(sqrt(p)-1)^2-1))*K2*Nz;

rhs_odd=term1-term2;

else

term1=1/(sqrt(p)+1)^t1*(aa1*(sqrt(p)+1)/(sqrt(p)-1)*sqrt(p)/(a*(sqrt(p)+1)^2-1))*K1*Nz;

term2=1/(sqrt(p)-1)^t1*(aa2*(sqrt(p)-1)/(sqrt(p)+1)*sqrt(p)/(a*(sqrt(p)-1)^2-1))*K2*Nz;

rhs_odd=term1+term2;

end

diff_odd=lhs_odd-rhs_odd;

check_error=checkmatrix(diff_odd);

if check_error>.5

disp('Error in Step 7, t odd.')

else

disp('Step 7 okay for odd t.')

end

%********************************************************************

%Case 2: t even

bb1=-psi1/(1-x1*psi1); %introduces abreviations bb1 and bb2.

bb2=-psi2/(1-x2*psi2);

%computes D_sigma(log p_t) from (13), t odd

% and computes D_sigma(log p_t^F) from (12) (drops constant terms from both)

a1_even=evl1^t2*J1*v1+evl2^t2*J2*v2;

a2_even=evl1^(t2-1)*P*v1+evl2^(t2-1)*P*v2+evl1^t2*M1*(E+evl1/a*D*C)*v1+evl2^t2*M2*(E+evl2/a*D*C)*v2;

lhs_even=a1_even-a2_even; %computes D_sigma(log p_t - log p_t^F), t odd

%computes Phi_1 and Phi_2 (these are produced in the first equation in the proof of Prop 3 (under Case 1).

%they are used in the followiung computation.

Phi1=(1+evl1^2/a)*(K1-L1)*v1-((1-evl1)/a+evl1^2/a^2)*L1*A*P*v1;

Phi2=(1+evl2^2/a)*(K2-L2)*v2-((1-evl2)/a+evl2^2/a^2)*L2*A*P*v2;

%computes D_sigma(log p_t - log p_t^F), t even, from cases (i) and (ii) in proof of Prop.3

if P_adj(z1+1,1)==1 & P_adj(z2+1,1)==1

yy1=1/(sqrt(p)+1)^(t2-1)*sqrt(p)*(p-1)/(2*(a*(sqrt(p)+1)^2-1))*(bb1*1/(a*(sqrt(p)+1)^2)-a)*F;

yy2=1/(sqrt(p)-1)^(t2-1)*sqrt(p)*(p-1)/(2*(a*(sqrt(p)-1)^2-1))*(bb2*1/(a*(sqrt(p)-1)^2)-a)*F;

yy3=1/a*A*(evl1^(t2+1)*Phi1+evl2^(t2+1)*Phi2);

rhs_even=yy1+yy2+yy3;

else

yy1=1/(sqrt(p)+1)^(t2-1)*sqrt(p)*(p-1)/(2*(a*(sqrt(p)+1)^2-1))*((sqrt(p)+1)/(sqrt(p)-1))*(-bb1*1/(a*(sqrt(p)+1)^2)+a)*F;

yy2=1/(sqrt(p)-1)^(t2-1)*sqrt(p)*(p-1)/(2*(a*(sqrt(p)-1)^2-1))*((sqrt(p)-1)/(sqrt(p)+1))*(-bb2*1/(a*(sqrt(p)-1)^2)+a)*F;

yy3=1/a*A*(evl1^(t2+1)*Phi1+evl2^(t2+1)*Phi2);

rhs_even=yy1+yy2+yy3;

end

diff_even=lhs_even-rhs_even;

check_error=checkmatrix(diff_even);

if check_error>.5

disp('Error in Step 7, t even.')

else

disp('Step 7 okay for even t.')

end

%**********************************************************************

%Step 8: Check psi^hat

%**********************************************************************

%Compute psi^hat according to def'n in Section 7.4.4

lhs_psiH1=mu1+eta1*((p-5)/4-1);

lhs_psiH2=mu2+eta2*((p-5)/4-1);

%Checking bounds on psi^hat.

if lhs_psiH1<=-1 | lhs_psiH1>=0

disp('Something wrong with psi1^hat')

end

if lhs_psiH2<=-1 | lhs_psiH2>=0

disp('Something wrong with psi2^hat')

end

%Compute psi_hat according to (35) (using polynomials of Section 7.4.4)

n1=a^4*(sqrt(p)+1)^4*(p-1)^2-a^3*(sqrt(p)+1)^4*(p-1)*(p-7)-a^2*(sqrt(p)+1)^2*(7*p-15)+a*(sqrt(p)+1)^2*(p-13)+4;

d1=a^4*(sqrt(p)+1)^4*(p-1)^2-2*a^2*(sqrt(p)+1)^2*(p+1)+1;

n2=a^4*(sqrt(p)-1)^4*(p-1)^2-a^3*(sqrt(p)-1)^4*(p-1)*(p-7)-a^2*(sqrt(p)-1)^2*(7*p-15)+a*(sqrt(p)-1)^2*(p-13)+4;

d2=a^4*(sqrt(p)-1)^4*(p-1)^2-2*a^2*(sqrt(p)-1)^2*(p+1)+1;

rhs_psiH1=-(a^2*(sqrt(p)+1)^2)/(4*(a^2*(sqrt(p)+1)^2-1))*n1/d1;

rhs_psiH2=-(a^2*(sqrt(p)-1)^2)/(4*(a^2*(sqrt(p)-1)^2-1))*n2/d2;

diff_psiH1=lhs_psiH1-rhs_psiH1;

diff_psiH2=lhs_psiH2-rhs_psiH2;

check_error1=checkmatrix(diff_psiH1);

check_error2=checkmatrix(diff_psiH2);

if check_error1+check_error2>.5

disp('Error in Step 8.')

else

disp('Step 8 okay.')

end