Replication of “A model of competitive stock trading volume”

I am replicating the paper "A model of competitive stock trading volume" by Jiang Wang http://web.mit.edu/wangj/www/pap/Wang94.pdf for a research project . I realize this is a very specific question but since this is not a very complicated paper and it is quite highly cited I am hoping maybe someone else has also replicated it and can help me out.

As far as I can tell solving the model boils down to

1. Solving a Kalman filtering problem (static Ricatti equation since we are interested in steady state solutions)
2. Obtaining the symmetric matrix of constants v in the optimal value function for the informed/uninformed investor (finding zeros of a system of equations)
3. Using the Kalman gain matrix and v from step 2 to find the values of prices p0, pF, pZ (finding zeros of another system of equations)

However, when I carry out steps 1 and 2 above for a given set of values p0, pF and pZ I cannot back out those values when performing step 3.

I have re-derived all the math to make sure there are no typos and there were none (in the solution process appendix) apart from the specification for the constant $$a_\theta$$ which is the coefficient of the AR(1) process describing the error that uninformed traders make in estimating the fundamental component of the stock dividend $$\hat{F}_t - F_t = a_{\theta}(\hat{F}_{t-1} - F_{t-1}) + b_{\theta} \epsilon_t$$. The Matlab script I am using (should run as is on Matlab 2016b or higher otherwise the eqPrices(v_u, v_i, k, p) function should be put in a separate file) is listed below:

% Parameters from Fig 1
r  = 0.2; gamma = 1.0; beta = 0.9;
aF = 0.9; aZ = 0.9;
varD = 2.0; varF = 2.0; varZ = 1.5; varQ = 1.0; covarDQ = 0.6;
omega = 0.1;

R = 1 + r;
a = aF/(R-aF);
% Parameters included to test solution alg (should be able to back
% these out when we solve for p0, pF, pZ with the optimal investment
% functions)
varS = 1.0;
p0 = 10;
pF = 2;
pZ = 4;

% vectors
b_D = [1, 0, 0, 0, 0]; b_F = [0, 1, 0, 0, 0]; b_Z = [0, 0, 1, 0, 0]; b_q = [0, 0, 0, 1, 0]; b_S = [0, 0, 0, 0, 1];

% state space and obs eqs
a_z = [aF 0; 0 aZ];
b_z = [b_F; b_D];
a_s = [1, 0; 1, 0; pF, pZ];
b_s = [b_S; b_D; zeros(1,5)];

% matrices
eps = [varD varF varZ varQ varS]; S = diag(eps);
S(1,4) = covarDQ; S(4,1) = S(1,4);
Szz = b_z*S*b_z';
Sss = b_s*S*b_s';
% Step 1
% Estimate stationary cond cov matrix O with dare() [dis alg Riccati eq]
% Check https://www.mathworks.com/help/control/ref/dare.html for reference
A_dare = a_z;
E_dare = eye(2);
B_dare = a_z'*a_s';
Q_dare = Szz;
R_dare = a_s*Szz*a_s'+Sss;
S_dare = Szz*a_s';
[O_dare,L_dare,G_dare] = dare(A_dare,B_dare,Q_dare,R_dare,S_dare,E_dare);
k = (a_z*O_dare*a_z' + Szz)*a_s'/(B_dare'*O_dare*B_dare+R_dare);
% Note that k = G_dare'

a_theta = aF - aF*k(1,1) - aF*k(1,2) - pF*k(1,3)*(aF+aZ);
b_theta = [k(1,2), k(1,1)+k(1,2)+k(1,3)*pF - 1, k(1,3)*pZ, 0, k(1,1)];

% Excess return coefficients
e0 = r*p0;
eZ = (R-aZ)*pZ;
e_theta = (R - a_theta)*(a-pF);
b_Q = (1+a)*b_F + b_D - pZ*b_Z + (a-pF)*b_theta;

%%%%%%%%%%%%%%%% Uninformed investor
% State vector variables
a_psi_u = [1 0;0 aZ];
b_psi_u = [zeros(1,6); -aZ*(pF/pZ) + a_theta*(pF/pZ), b_Z + b_theta*(pF/pZ)];
% Excess return vector variables
e_2_u = [e0, eZ];
b_2_u = [-e_theta - (pF/pZ)*eZ, b_Q];
% Value function variables
v_aa_u = @(v) a_psi_u'*[v(1) v(2); v(2) v(3)]*a_psi_u;
v_bb_u = @(v) b_psi_u'*[v(1) v(2); v(2) v(3)]*b_psi_u;
v_ab_u = @(v) a_psi_u'*[v(1) v(2); v(2) v(3)]*b_psi_u;
S_u = [b_theta*S*b_theta' zeros(1,5); zeros(5,1) S];
Omega_u = @(v) inv(inv(S_u)+v_bb_u(v));
Gamma_u = @(v) inv(b_2_u*Omega_u(v)*b_2_u');
g_u = @(v) e_2_u - b_2_u*Omega_u(v)*v_ab_u(v)';
d_u = @(v) det(inv(Omega_u(v))*S_u)^(-1/2);
m_u = @(v) v_aa_u(v) - v_ab_u(v)*Omega_u(v)*v_ab_u(v)' + g_u(v)'*Gamma_u(v)*g_u(v);
% Variables related to equation for constant term v
index_mat_u = zeros(2,2);
index_mat_u(1,1) = 1;
c_bar_u = @(v) -1/(gamma*R)*log(r*beta*d_u(v));

%%%%%%%%%%%%%%%% Informed investor
% State vector variables
a_psi_i = [1 0 0; 0 aZ 0; 0 0 a_theta];
b_psi_i = [zeros(1,5); b_Z; b_theta];
% Excess return vector variables
e_2_i = [e0 eZ e_theta; 0 1 0];
b_2_i = [b_Q; b_q];
% Value function variables
v_aa_i = @(v) a_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*a_psi_i;
v_bb_i = @(v) b_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*b_psi_i;
v_ab_i = @(v) a_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*b_psi_i;
S_i = S;
Omega_i = @(v) inv(inv(S_i)+v_bb_i(v));
Gamma_i = @(v) inv(b_2_i*Omega_i(v)*b_2_i');
g_i = @(v) e_2_i - b_2_i*Omega_i(v)*v_ab_i(v)';
d_i = @(v) det(inv(Omega_i(v))*S_i)^(-1/2);
m_i = @(v) v_aa_i(v) - v_ab_i(v)*Omega_i(v)*v_ab_i(v)' + g_i(v)'*Gamma_i(v)*g_i(v);
% Variables related to equation for constant term v
index_mat_i = zeros(3,3);
index_mat_i(1,1) = 1;
c_bar_i = @(v) -1/(gamma*R)*log(r*beta*d_i(v));

% Step 2
% Solve for matrices of constants v_i and v_u (from value function ansatz)
v_u = fsolve(@(v) (1/R)*m_u(v) - [v(1) v(2); v(2) v(3)] + (gamma*c_bar_u(v)+log(r/R))*index_mat_u, [0,0,0]);
v_i = fsolve(@(v) (1/R)*m_i(v) - [v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)] + (gamma*c_bar_i(v) + log(r/R))*index_mat_i, [0,0,0,0,0,0]);

% Step 3
% After obtaining k from the filtering problem and v from the opt problem
% we can solve for p0, pF, pZ
p = fsolve(@(p) eqPrices(v_u, v_i, k, p), [1;1;1])

function F = eqPrices(v_u, v_i, k, p)
% Parameters from Fig 1
r  = 0.2; gamma = 1.0; beta = 0.9;
aF = 0.9; aZ = 0.9;
varD = 2.0; varF = 2.0; varZ = 1.5; varQ = 1.0; covarDQ = 0.6;
omega = 0.1;

R = 1 + r;
a = aF/(R-aF);
% Same varS as in Step 1 and Step 2
varS = 1.0;

% vectors
b_D = [1, 0, 0, 0, 0]; b_F = [0, 1, 0, 0, 0]; b_Z = [0, 0, 1, 0, 0]; b_q = [0, 0, 0, 1, 0]; b_S = [0, 0, 0, 0, 1];

% state space and obs eqs
a_z = [aF 0; 0 aZ];
b_z = [b_F; b_D];
a_s = [1, 0; 1, 0; p(2), p(3)];
b_s = [b_S; b_D; zeros(1,5)];

% matrices
eps = [varD varF varZ varQ varS]; S = diag(eps);
S(1,4) = covarDQ; S(4,1) = S(1,4);
Szz = b_z*S*b_z';
Sss = b_s*S*b_s';

a_theta = aF - aF*k(1,1) - aF*k(1,2) - p(2)*k(1,3)*(aF+aZ);
b_theta = [k(1,2), k(1,1)+k(1,2)+k(1,3)*p(2)-1, k(1,3)*p(3), 0, k(1,1)];

% Excess return coefficients
e0 = r*p(1);
eZ = (R-aZ)*p(3);
e_theta = (R - a_theta)*(a-p(2));
b_Q = (1+a)*b_F + b_D - p(3)*b_Z + (a-p(2))*b_theta;

%%%%%%%%%%%%%%%% Uninformed investor
% State vector variables
a_psi_u = [1 0;0 aZ];
b_psi_u = [zeros(1,6); -aZ*(p(2)/p(3)) + a_theta*(p(2)/p(3)), b_Z + (p(2)/p(3))*b_theta];
% Excess return vector variables
e_2_u = [e0, eZ];
b_2_u = [-e_theta - (p(2)/p(3))*eZ, b_Q];
% Value function variables
v_aa_u = @(v) a_psi_u'*[v(1) v(2); v(2) v(3)]*a_psi_u;
v_bb_u = @(v) b_psi_u'*[v(1) v(2); v(2) v(3)]*b_psi_u;
v_ab_u = @(v) a_psi_u'*[v(1) v(2); v(2) v(3)]*b_psi_u;
S_u = [b_theta*S*b_theta' zeros(1,5); zeros(5,1) S];
Omega_u = @(v) inv(inv(S_u)+v_bb_u(v));
Gamma_u = @(v) inv(b_2_u*Omega_u(v)*b_2_u');
g_u = @(v) e_2_u - b_2_u*Omega_u(v)*v_ab_u(v)';
d_u = @(v) det(inv(Omega_u(v))*S_u)^(-1/2);
m_u = @(v) v_aa_u(v) - v_ab_u(v)*Omega_u(v)*v_ab_u(v)' + g_u(v)'*Gamma_u(v)*g_u(v);
% Variables related to equation for constant term v
index_mat_u = zeros(2,2);
index_mat_u(1,1) = 1;
c_bar_u = @(v) -1/(gamma*R)*log(r*beta*d_u(v));

%%%%%%%%%%%%%%%% Informed investor
% State vector variables
a_psi_i = [1 0 0; 0 aZ 0; 0 0 a_theta];
b_psi_i = [zeros(1,5); b_Z; b_theta];
% Excess return vector variables
e_2_i = [e0 eZ e_theta; 0 1 0];
b_2_i = [b_Q; b_q];
% Value function variables
v_aa_i = @(v) a_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*a_psi_i;
v_bb_i = @(v) b_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*b_psi_i;
v_ab_i = @(v) a_psi_i'*[v(1) v(2) v(3); v(2) v(4) v(5); v(3) v(5) v(6)]*b_psi_i;
S_i = S;
Omega_i = @(v) inv(inv(S_i)+v_bb_i(v));
Gamma_i = @(v) inv(b_2_i*Omega_i(v)*b_2_i');
g_i = @(v) e_2_i - b_2_i*Omega_i(v)*v_ab_i(v)';
d_i = @(v) det(inv(Omega_i(v))*S_i)^(-1/2);
m_i = @(v) v_aa_i(v) - v_ab_i(v)*Omega_i(v)*v_ab_i(v)' + g_i(v)'*Gamma_i(v)*g_i(v);
% Variables related to equation for constant term v
index_mat_i = zeros(3,3);
index_mat_i(1,1) = 1;
c_bar_i = @(v) -1/(gamma*R)*log(r*beta*d_i(v));

alpha = r*gamma/R;
X_t_u = (1/alpha)*Gamma_u(v_u)*(e_2_u - b_2_u*Omega_u(v_u)*v_ab_u(v_u)');
X_t_i = (1/alpha)*Gamma_i(v_i)*(e_2_i - b_2_i*Omega_i(v_i)*v_ab_i(v_i)');

F = [omega*X_t_i(1,1) + (1-omega)*X_t_u(1,1) - 1;
omega*X_t_i(1,2) + (1-omega)*X_t_u(1,2);
omega*X_t_i(1,3) + (1-omega)*(p(2)/p(3))*X_t_u(1,2)];
end

Where the function eqPrices(v_u, v_i, k, p) takes in k, v_u, v_i from Step 1 and 2 as given and minimizes for the market clearing condition which translates to finding the zeros of:

F = [omega*X_t_i(1,1) + (1-omega)*X_t_u(1,1) - 1;
omega*X_t_i(1,2) + (1-omega)*X_t_u(1,2);
omega*X_t_i(1,3) + (1-omega)*(p(2)/p(3))*X_t_u(1,2)];

Here X_t_i and X_t_u are the optimal investment functions for the informed and uninformed investor respectively (a 1x3 and a 1x2 vector of constants). Omega is the share of informed investors.

The only tricky part for me was to recognize that the vector of state variables for the uninformed investor can be written as:

where $$\Theta_{t-1} =\hat{F}_{t-1} - F_{t-1}$$, $$a_{\theta}$$ and $$b_{\theta}$$ are functions of the steady state Kalman gain matrix coefficients and p0, pF and pZ. Their expressions are taken from the code above and are given as:

a_theta = aF - aF*k(1,1) - aF*k(1,2) - pF*k(1,3)*(aF+aZ);
b_theta = [k(1,2), k(1,1)+k(1,2)+k(1,3)*pF - 1, k(1,3)*pZ, 0, k(1,1)];

Any suggestions on why I can't back out p0, pF and pZ or on any fundamental mistakes in my understanding of the solution process are highly appreciated.