vpasolve returns Empty sym: 0-by-1

I have the following code:

%**************************************************%
%       Parameters of the cost function            %
%**************************************************%
A = 5;
alpha = 0.5;
beta = 0.5;
gama = 0.6;
%**************************************************%
%       Distribution of types                      %
%***************************************************
pll = 0.36;
phl = 0.14;
plh = 0.14;
phh = 0.36;
%**************************************************%
%      Asymmetric information parameters           %
%**************************************************%
theta_qL = 2;
theta_qH = 4;
theta_tL= 4;
theta_tH= 8;
%**************************************************%
%        Deltas                                    %
%***************************************************
delta_HH_LH = ((theta_qH)^beta)*((theta_tH)^(1-beta))...
    -((theta_qL)^beta)*((theta_tH)^(1-beta));
delta_HH_HL = ((theta_qH)^beta)*((theta_tH)^(1-beta))...
    -((theta_qH)^beta)*((theta_tL)^(1-beta));
delta_HH_LL = ((theta_qH)^beta)*((theta_tH)^(1-beta))...
    -((theta_qL)^beta)*((theta_tL)^(1-beta));
delta_LH_LL = ((theta_qL)^beta)*((theta_tH)^(1-beta))...
    -((theta_qL)^beta)*((theta_tL)^(1-beta));
delta_HL_LL = ((theta_qH)^beta)*((theta_tL)^(1-beta))...
    -((theta_qL)^beta)*((theta_tL)^(1-beta));
%**************************************************%

K1 = gama*(plh*delta_HH_LH + phl*delta_HH_HL + pll* delta_HH_LL);
K2 = gama*((phl+plh)*(delta_HH_LL)- phl*delta_HH_HL - plh*delta_HH_LH);
K3 = gama*((theta_qH)^beta)*((theta_tH)^(1-beta));
K4 = gama*delta_HH_LH;
K5 = gama*delta_HH_HL;
K6 = gama*delta_HH_LL;

Q1 = gama*((theta_qL)^beta)*((theta_tH)^(1-beta));
Q2 = gama*delta_LH_LL;

R1 = gama*((theta_qH)^beta)*((theta_tL)^(1-beta));
R2 = gama*delta_HL_LL;

S1 = gama*((theta_qL)^beta)*((theta_tL)^(1-beta));

%**************************************************%
%                Computing qhh                     %
%**************************************************%

thh = 2;

syms qhh

vpasolve(-(A*alpha*(qhh^(alpha - 1))*(thh^(alpha - 1)))...
       + (((qhh*thh)^(gama - 1))*thh)* K3 ...
       + (((((qhh*thh)^(gama - 1))*thh)*K1/((phh/plh)*(phh + (((qhh*thh)^(gama - 1))*thh)*K2))...
       + (plh/phh))*((((qhh*thh)^(gama - 1))*thh)*K4)) + (((((qhh*thh)^(gama - 1))*thh)*K1/((phh/phl)*(phh + (((qhh*thh)^(gama - 1))*thh)*K2))...
       + (phl/phh))*((((qhh*thh)^(gama - 1))*thh)*K5)) + ((pll/phh) - ((((qhh*thh)^(gama - 1))*thh)* K1/(phh^2 + (((qhh*thh)^(gama - 1))*thh)*K2)))*((((qhh*thh)^(gama - 1))*thh)*K6) == 0,qhh);

As a result I obtain the "satisfactory" result:

ans =

0.034331822356251028338781025191365

However, as soon as I change the values in the distribution of types to anything below 0.36 for pll and phh, and 0.14 for plh and phl, I obtain Empty sym: 0-by-1.

For example, in the code above, use:

%**************************************************%
%       Distribution of types                      %
%***************************************************
pll = 0.35;
phl = 0.15;
plh = 0.15;
phh = 0.35;

And I get:

ans =

Empty sym: 0-by-1

Question: What is the reason for this? My model, of course, is very sensitive to the distribution of types but I would expect vpasolve to return a value, not to tell me it cannot find a solution.

Thanks for your help,

Luis

[John D'Errico]

And, do you know, positively, that there is a solution for that changed set of parameters? Perhaps the fact that no solution was found might be a suggestion that no solution exists?

[Luis Boscan]

Hi John, according to my model yes. However, the solution with pll = 0.36; phl = 0.14; plh = 0.14; phh = 0.36; is close to the theoretical properties I expect. That is, there should be a set of parameters for which the answer is even closer to zero.

Answer 1

[John D'Errico]

Well, you did not answer my question about KNOWING a solution exists. Close enough to zero is not zero. Anyway, this is a ONE variable problem.

PLOT EVERYTHING. If you can plot it, before you just throw a solver at something, plot it.

Eqn = -(A*alpha*(qhh^(alpha - 1))*(thh^(alpha - 1)))...
     + (((qhh*thh)^(gama - 1))*thh)* K3 ...
     + (((((qhh*thh)^(gama - 1))*thh)*K1/((phh/plh)*(phh + (((qhh*thh)^(gama - 1))*thh)*K2))...
     + (plh/phh))*((((qhh*thh)^(gama - 1))*thh)*K4)) + (((((qhh*thh)^(gama - 1))*thh)*K1/((phh/phl)*(phh + (((qhh*thh)^(gama - 1))*thh)*K2))...
     + (phl/phh))*((((qhh*thh)^(gama - 1))*thh)*K5)) + ((pll/phh) - ((((qhh*thh)^(gama - 1))*thh)* K1/(phh^2 + (((qhh*thh)^(gama - 1))*thh)*K2)))*((((qhh*thh)^(gama - 1))*thh)*K6);

So, knowing that you expect the solution to be near zero, try this plot:

ezplot(Eqn,[0,.1])

Clearly the upper branch of this will diverge, never crossing zero. So zoom in on the lower branch.

ezplot(Eqn,[0,0.01])
grid on

So, lets see if vpasolve is getting lost, searching in that divergent domain.

vpasolve(Eqn,0.01)
ans =
Empty sym: 0-by-1

vpasolve(Eqn,0.001)
ans =
0.0011108227632165429244076492686323

1. Just because you want a solution to exist does not mean it will exist. It MIGHT exist, but sometimes you need to do a little extra thinking to make things work.

2. When all else fails, MAKE A PLOT. In fact, even before you throw a messy problem at a solver, MAKE A PLOT. Plot everything that you can plot. You can often learn something from those plots.