How to avoid NaN with numerical integration (quadrature)

Hi, I have written a program that I divided in several functions.

"*puzzle*" is the main function. But my question is related to a problem I have with quadv(Look in "*LnpdFut*" function).

s_ is a parameter of "*pdFut*" that I passed as a global to "*Lnpdfut*" function, the argument of integration being v.

In the loop in "*LnpdFut*", s_ takes successively values of the vector sglob_, but for the first two values of sglob_, it doesn't work. I get NaN.

Could someone help me figure out the reason? I know it's long but just run it, you can see for yourself.

    function puzzle

global lnpdc_ sglob_ gamma_ g_ G_ rf_ phi_ sigmav_ sigmaw_ rho_ Sbar_ sbar_ smax_ delta_ LimInf_ LimSup_ ndivs_

gamma_= 2;                                       
g_= 0.0174;                                     
G_= exp(g_);
rf_ = 0.023;                                          
phi_= 0.9016;                                    
sigmav_= 0.01186;                                
sigmaw_= 0.465;                                  
rho_= 0.3;
Sbar_= sigmav_*(gamma_/(1-phi_))^(0.5);            
sbar_= log(Sbar_);
smax_= sbar_ + (1-Sbar_^2)/2;                     
delta_ = 1/exp(rf_)*G_^(gamma_)*exp(-gamma_*(1-phi_)/2); 
ndivs_ = 15;
LimInf_ = -3*sigmav_;
LimSup_ = 3*sigmav_;

%Main function puzzle body

sglob_ = Setgrid(smax_,sbar_,ndivs_); % Grille
 fprintf('sglob_ est %d x %d\n', size(sglob_,1),size(sglob_,2));   
lnpdc_ = zeros(length(sglob_),1);
pdcInt = Lnpdfut(sglob_); % intègre numériquement
fprintf('pdcInt:%d x %d\n',size(pdcInt,1),size(pdcInt,2));
disp(pdcInt);
%==================================================================

function k = Lnpdfut(sg)
%Calls Futpd (below) and uses quadv

global s_ nwpdc_

nwpdc_ = zeros(length(lnpdc_),1);

j=1;
    while(j<=length(lnpdc_))
    s_= sg(j);                         
nwpdc_(j)= quadv(@Futpd,LimInf_,LimSup_); %integration sur v

j=j+1;
    end
k = nwpdc_; 

end

%------------------------------------------------------------------

function pdFut = Futpd(v)
%calls Transit belows
global s_

pdFut = delta_*G_^(1-gamma_).*exp(-gamma_*(Transit(s_,v)-s_.*ones(length(Transit(s_,v)),1))).*...  
    exp((1-gamma_)*v).*(1 +     exp(interp1q(sglob_,lnpdc_,Transit(s_,v)))).*pdf('norm',v,0,sigmav_);
end

%------------------------------------------------------------------

function Futstate = Transit(s,v) 
%Calls Lam below

Futstate =(1-phi_)*sbar_ + phi_.*s + Lam(s).*v;

if max(Futstate >=smax_)==1,
    Futstate = (Futstate<smax_).*Futstate + (Futstate>=smax_).*smax_; 
end
end
%------------------------------------------------------------------

function lambda = Lam(s) 

express= 1-2*(s-sbar_);
express =(express>=0).*express;  
express= express.^(0.5);
express= (1./Sbar_).*express-1;
express= (s<=smax_).*express; 

lambda = express;
end

%------------------------------------------------------------------

function s_etat = Setgrid(smx,sb,nd)

s_etat = 0:exp(smx)/nd:exp(smx); %donne (ndivs+1) points
s_etat = log(s_etat(2:end));
if max(smx==s_etat)== 0,  
 s_etat = horzcat(s_etat,smx);  
 end

if max(sb==s_etat)== 0,
  s_etat = horzcat(s_etat,sb);  
end

s_etat = sort(s_etat);
s_etat = s_etat';

end
end

[Matt Fig]

I would try it, except you forgot to use the {} Code button when pasting your code so that it looks like a mess. You can go back and edit the question so that the rest of us can simply past it into MATLAB.

[Ismael]

It's done. Thanks for pointing this to me.

I'm looking forward to your suggestion(s)

Answer 1

[Andrew Newell]

If you set a debugging breakpoint at the call where the warning occurs, then

Futpd(LimInf_)

ans = 

     NaN

EDIT: The NaN occurs in this expression with v=LimInf_:

interp1q(sglob_,lnpdc_,Transit(s_,v))

The problem is that you are trying to extrapolate data for x between about -5 and -2 to x=-223, and interp1q doesn't like that. If you are always extrapolating, you could try

interp1(sglob_,lnpdc_,Transit(s_,v),'linear','extrap')

This would work for this case, where lnpdc_ is zero for all values of x, but in general extrapolation is dangerous.

Programming note: Don't use global variables.

Answer 2

[bym]

I have no idea, but don't you have to declare the global variables in each subfunction? So for

Lnpdfut(sg)

you would also have to declare

LimInf_,LimSup_

[Ismael]

Thank you!!! Your answer inspired me. I looked closely at interp in matlab help and I found an alternative solution to extrapolate using ppval(). I needed this for my thesis.

Thank you!!!