solve numerically a system of three equations the first of which defines a function and the rest of two involve integrals

I want to solve numerically the equilibrium characterised by the following system of three equations:

where are scalar parameters, Φ the standard normal CDF and the bivariate normal PDF. The goal is to compute numerical solutions by solving the above system of equations simultaneously for the function and the scalars .

I have been searching on the forum. When involving integral equations in the system, integral2() and fsolve() are recommended. I looked over the documentations, but cannot figure out how to include the first equation of the above system in the formulation.

Could anyone suggest a way out? I appreicate.

[Walter Roberson]

Could you confirm that theta is an unknown function and you want your numeric system to find that function? Because that is something that numeric systems cannot do. If theta is a function of known form with a couple of parameters to be found, that might possibly be workable.

[Yuequan]

Thanks for a prompt response. θ is a function of , but the form is unknown. I do not need the numeric system to find the (exact) function of θ. My interest concerns values of , but I guess I have to somehow "compute" in order to solve for ? In other words, a list of that satisfies the first equation (condition on ) should suffice.

PS. θ is strictly decreasing in and .

[David Goodmanson]

hello YG,

It appears that the domain of theta(zeta1,zeta2) is all of R^2. Could you elaborate on what 'strictly decreasing' means here?

[Walter Roberson]

Φ the standard normal CDF

with respect to which variable in each case? Is Φ being used as a multiplier, or is Φ(expression) applying CDF to the expression? CDF is something applied to random variables, but I do not see any random variables?

[Yuequan]

I mean the partial derivative is strictly negative, i.e. .

In addition, let and . Then, the partial derivative is also strictly negative, i.e. .

[Yuequan]

@Walter, Φ(.) is evaluated at and , respectively, in eqn 2 and 3. The original problem involves another standard normal random variable where Φ comes in.

[Walter Roberson]

I need to know which variable name to put into the call to the symbolic cumulative distribution function as the random variable being acted on. Which of the variables η or or or is the random variable?

[Yuequan]

@Walter, thanks. In eqn2 and eqn3, random variables are , which jointly follows a bivariate normal distribution. In eqn1, there is no random variables (i.e. no probability operation involved) although θ is a function of .

[David Goodmanson]

Hi YG, do you have suggested values for the parameters?

[Walter Roberson]

Equation 1 involves two cdf operations, so there are random variables involved

[David Goodmanson]

also, how complicated is the bivariate normal pdf? Is it as simple as

normalization_factor * exp(-(zeta1^2+zeta2^2)/2)   ?

[Yuequan]

@Walter, I don't know much about symbolic cumulative distribution function in Matlab so I might misunderstand what you mean by random variables. As far as I know, characterises the distributoin of a (standard normal) random variable, but is itself a deterministic function (not closed form though). So that's why I said eqn1 doesn't involve extra random variables. But if needed, I guess we can always write where ℙ is the probability operator and Z a random variable following a normal distribution with 0 mean and 1 variance.

Answer 1

[David Goodmanson]

Hello YG,

Fortunately the first equation is amenable to interation, at least for some parameters. In the code below,

zeta1 --> x,    zeta2 -->  y,    x1 --.> x0,    x2 --> y0.

You didn't say if the bivariate normal is general, but that's incorporated here. RIght now the paramters for that are contained in the function definition. Also inlcuded is separate constants C and D for x and y.

Clearly C,D <1 is necessary. The code has problems if sqrt(eta) is more than just a smidge above 1, but does all right for smaller values.

If you change delx and dely to .25 for plotting purposes, and something like

a = 1/3;  C = .8;  D = .8;  sn = 1;

with,say, a plain vanilla bivariate normal,

mux = 0;   muy = 0;   sigx = 1;  sigy = 1;   rho = 0;

and then do the commented-out surf(x,y), you can see a nice plot of the theta function dropping monotonically from 1 at (-6,6) to 0 at (6,6)

a = 4/5;
C = .88;
D = .5;
sn = 1;    % sn = sqrt(eta)

delx = .012;
dely = .012;
delxy = delx*dely;
xx = -6:delx:6;
yy = -6:dely:6;
[x y] = meshgrid(xx,yy);
f = bivnorm(x,y);       % parameters are contained within the function
% initial guesses
th = normcdf(x).*normcdf(y);
x0 = 0;
y0 = 0;
% tolerances
tolTH = 1e-5;
tolC = 1e-5;
tolD = 1e-5;


thold = th+2*tolTH;        % start off while loop
while  max(max(abs(thold-th))) > tolTH
 x0 = Gx(x,y,delxy,th,sn,x0,f,C,tolC);
 y0 = Gy(x,y,delxy,th,sn,y0,f,D,tolD);
 thold = th;
 % iterate theta
 th = a*normcdf(sn*((x0-x)-th)) +(1-a)*normcdf(sn*((y0-y)-th));
end

% verify, these should be small
Ccheck = trapz(trapz(normcdf(sn*(th-(x0-x))).*f))*delx*dely - C
Dcheck = trapz(trapz(normcdf(sn*(th-(y0-y))).*f))*delx*dely - D
thnew = a*normcdf(sn*((x0-x)-th)) +(1-a)*normcdf(sn*((y0-y)-th));
thcheck = max(max(abs((thnew-th))))
x0
y0
% surf(x,y,th)

function x0new = Gx(x,y,delxy,th,sn,x0,f,C,tolC)
% root for the second eqn
x0old = x0 +2*tolC;      % start off while loop
while abs(x0old-x0)>tolC
  A = trapz(trapz(normcdf(sn*(th-(x0-x))).*f))*delxy - C;
  derivA = -trapz(trapz(normpdf(sn*(th-(x0-x))).*f))*delxy;
  x0old = x0;
  x0 = x0 - A/derivA;
end
x0new = x0;
end

function y0new = Gy(x,y,delxy,th,sn,y0,f,D,tolD)
% root for the third eqn
y0old = y0 +2*tolD;     % start off while loop
while abs(y0old-y0)>tolD
  A = trapz(trapz(normcdf(sn*(th-(y0-y))).*f))*delxy - D;
  derivA = -trapz(trapz(normpdf(sn*(th-(y0-y))).*f))*delxy;
  y0old = y0;
  y0 = y0 - A/derivA; 
end
y0new = y0;
end

function f = bivnorm(x,y)
% bivariate normal pdf for meshgrid-style x and y matrices
mux = .1;
muy = -.1;
sigx = 1.2;
sigy = .88;
rho = .2;
xs = x - mux;
ys = y - muy;
rhoc = sqrt(1-rho^2);
A = xs.^2/sigx^2 -2*rho*xs.*ys/(sigx*sigy) + ys.^2/sigy^2;
f = 1/(2*pi*sigx*sigy*rhoc)*exp(-A/(2*rhoc^2));
end

,

[Yuequan]

@David, thanks! Sorry for a late response - I was working on the problem and didn't go on the forum for a few days. I wanted the bivariate normal density to be general. I am inexperinced with Matlab and I'm studying your codes. Here are a few questions:

1. why use "th = normcdf(x).*normcdf(y)" as the initial guess for theta?
2. in the local function Gx and Gy, the derivative of the function is also calculated. I guess it's for the choice of the new x and y, but why the form of "x0 = x0 - A/derivA"?

Also, with a lot of help from a friend, I came up with the following code, which makes use of fsolve(). The problem is that it takes 5 minutes to run for a (zeta1, zeta2) grid where each runs from -5 to 5 in an increment of 0.1. As I need to run the codes for a large number of parameter combinations, any suggestions about how to speed it up are appreciated. Thanks for your time!

% structural parameters
alpha = 0.5;
oppc = 0.5;
gamma1 = 0.5;
gamma2 = 0.5;
eta = 0.7;


% algorithm #3
% numerical parameters
z = 5;
z1min = -z;
z1max = z;
z2min = -z;
z2max = z;
dz = 0.1;

[z1, z2] = meshgrid(z1min:dz:z1max, z2min:dz:z2max);
g_info = [z1(:) z2(:)]; Nz = length(g_info);
zeta1 = g_info(:,1); zeta2 = g_info(:,2);


% parameters of interests
Rho = -1:0.05:1;
writematrix([],"equilibrium_solution.csv");

for rho = Rho
	%% iterative paramters
	iter = 0;
	max_iter = 100;
	relax = 0.02;

	error0 = 100; error1 = 100; error2 = 100;
	tolerance = 10e-3;


	%% initial guess
	t0 = repelem(0.5,Nz);
	x10 = 0.5;
	x20 = 0.5;

	options = optimoptions("fsolve","OptimalityTolerance",1e-10,"FunctionTolerance",1e-10,"StepTolerance",1e-10,"Display",'final-detailed'); %,'PlotFcn',@optimplotfirstorderopt, 'Algorithm','levenberg-marquardt'

	while iter < max_iter && (error0 > tolerance || error1 > tolerance || error2 > tolerance)
		%%% solve the equilibrium
		disp(['iteraton = ', num2str(iter)]);
		f = @(X)equilibrium(X, zeta1, zeta2, alpha, g_info, dz, oppc, gamma1, gamma2, rho, eta, Nz);
		X0 = [x10 x20 t0];
		[X1, f_diff1] = fsolve(f, X0, options); %
		f_diff0 = f(X0);

		%%% checking whether the solution converges
 		error0 = max(abs(X1(3:Nz+2) - X0(3:Nz+2))); error1 = max(abs(X1(1) - X0(1))); error2 = max(abs(X1(2)-X0(2)));
 		disp(['error0 = ', num2str(error0), '; error1 = ', num2str(error1), '; error2 = ', num2str(error2)]);
 
 		X_current = [X0(1) X0(2) X0(3:Nz+2)];
 
 		t0 = (1-relax).*X1(3:Nz+2) + relax.*X0(3:Nz+2);
 		x10 = (1-relax)*X1(1) + relax*X0(1);
 		x20 = (1-relax)*X1(2) + relax*X0(2);

		iter = iter+1;
	end

	%% check why the while loop ends
	if iter < max_iter && error0 < tolerance && error1 < tolerance && error2 < tolerance
	    disp('congratulations')
	else
	    disp('sorry')
    end

	%% assemble the solutions corresponding to parameters of interests
	writematrix(X_current, 'equilibrium_solution.csv', 'WriteMode', 'append');

end

toc

% plots
m = readmatrix("equilibrium_solution.csv");
m = [transpose(Rho),m];
plot(m(:,1), m(:,2));
xline(0);
axis([-1 1 0.475 0.5]);
saveas(gcf,"equilibrium_solution.png");

% local functions
function F = equilibrium(X, zeta1, zeta2, alpha, g_info, dz, oppc, gamma1, gamma2, rho, eta, Nz)

	% solve for unconditional 
	weight = mvnpdf(g_info, [], [(1/gamma1) (rho * sqrt(1/(gamma1 * gamma2))); (rho * sqrt(1/(gamma1 * gamma2))) (1/gamma2)]) * (dz * dz);
	v1 = normcdf(sqrt(eta) .* (transpose(X(3:Nz+2)) - X(1) + zeta1));
	v2 = normcdf(sqrt(eta) .* (transpose(X(3:Nz+2)) - X(2) + zeta2));

	F(1) = transpose(v1) * weight - oppc;
	F(2) = transpose(v2) * weight - oppc;

	% solve for conditional 
	F(3:Nz+2) = alpha .* normcdf(sqrt(eta) .* (transpose(X(3:Nz+2)) - X(1) + zeta1)) + (1-alpha) .* normcdf(sqrt(eta) .* (transpose(X(3:Nz+2)) - X(2) + zeta2)) - (1-transpose(X(3:Nz+2)));

end

Answer 2

[David Goodmanson]

Hi YG,

Not having the toolbox for fsolve, I can't do any comparison calculations. As is usually the case, there is a tradoff betweeen speed and precision. What kind of precision are you looking for? Compared to the method I used (and there could be other methods as well), I don't think that fsolve is necessarily better. It's clearly time consuming, and that is with a tolerance of 1e-2 thrown in which is pretty high. As I said I can't do comparisons, but I assume that your 5 min time is for the 41 instances of Rho, with a 101x101 grid, is that correct? The code I wrote repeated 41 times with 241x241 grid takes 5 sec, so there is a speed advantage if the precision can be shown to be at least as good. It's an interesting problem and I am going to keep going on it for awhile.

Regardless of the method, when x1 or x2 become large, then the grid in (zeta1, zeta2) is going to have problems because the function is shifted and may not fall to small enough values at the edge of the grid. Again that is going to depend on what kind of precision is wanted.

[1] I used th = normcdf(x).*normcdf(y) because once I realized that theta is really a cdf, I just wanted to start with a 2d cdf of some kind. Except that theta is a compementary cdf, 1 --> 0, not the usual cdf. I should have used th = normcdf(-x).*normcdf(-y) but it didn't make any difference. The iteration seems to be extremely robust to starting conditions, so now for an nxn grid I am just using ones(n,n).

[2] x0 = x0 - A/derivA is just Newton's method for root finding. I like the method because you don't need an upper and lower estimate for the root and if you have a good set of starting points you can calculate a whole lot of roots at once, such as the first 50 roots of a the bessel function. The method does require the derivative. Theta is numeric so usually the method would not apply, but this is an unusual situation where the derivative is also available at each point, numerically to the same level of accuracy as the function itself.