I am trying to solve an equation, but I am not getting the correct value so was wondering if anyone could see my mistake.

I am trying to solve two equations(Please read the attatched file for a better representation of what I am saying below):

1. The first, which I am trying to solve for theta1 is, g(theta1) = -log(kappa)/lambda, where g(s) = 2 [ Φ(a) + (s − 1)Φ(a(1 − s)) + [exp(− a^2 /2 − (a(1 − s))^2 / 2 ]/(a√ 2π) ]− s and Φ = phi is the cumalitive normal distribution function.
2. The second, which I am trying to solve for theta2 is, 1 = theta2 + (1 - kappa)/(kappa lambda h(theta2)), where h(t) = [2Φ(sigma(1-t)) - 1] and again Φ = phi is the cumalitive normal distribution function.

I have tried the following code:

beta = 1000

beta = 1000

alpha = 467

alpha = 467

kappa = 1 - alpha/beta % Influences the gamma* graphs, not sure how yet look at.

kappa = 0.5330

sigma = 4 % Appears to change the gamma graphs gradient/where it is at t = 1. It

sigma = 4

% also changes the value of theta1 (first intersection). The
% larger the value of sigma and thus the larger the value of a,
% the closer in value theta1 is to -log(kepa)/lambda.
a = sigma/2

a = 2

lambda = 3.2 % Appears to change the gamma graphs gradient/where it is at t = 1. It

lambda = 3.2000

% also changes the value of theta1 (first intersection).
b = -log(kappa)/lambda

b = 0.1966

Calculating θ_1 :

syms theta1
gtheta1 = 2*(normcdf(a) + (theta1 - 1) * normcdf(a*(1-theta1)) + 1/(a*sqrt(2*pi)) * (exp(-(a.^2)/2) - ...
    exp(-(a.^2*(1-theta1).^2)/2))) - theta1 == b;
soltheta1 = vpasolve(gtheta1, theta1)

soltheta1 = 

0.21283001670683141454156213680948

w = round(soltheta1, 3);
w1 = w + 0.001;
w2 = 1000 * w1;
w3 = (w + 0.002) * 1000;

Calculating θ_2:

syms theta2
gtheta2 = theta2 + (1 - kappa)/(kappa * lambda * (2*normcdf(sigma*(1 - theta2) / 2 ) - 1)) == 1;
soltheta2 = vpasolve(gtheta2, theta2)

soltheta2 = 

0.55954284069178087765059472947352

w = round(soltheta2, 3);
x = w + 0.001;
y = 1000 * x;
z = (w + 0.002) * 1000;

Note that the answers should be, theta1 = 0.175 and theta2 = 0.85

[Paul]

I got the same answer for theta_1. Why do you think it's incorrect?

[Charles]

This is because this should be the point of stationarity for a graph but it is not, this occurs at theta1 = 0.175. My professor also computes it to be that.

[Paul]

Perhaps double check the equation in Q1 of the attached pdf and/or the constants beta, alpha, and sigrma. A plot shows two roots, one at ~0.2128 and another at ~1.785

[Charles]

I think the equation must be wrong but I don't know how, the parameters values are definitely correct. The odd thing is that it works in another scenario to generate a correct theta1.