Double integral with complex numbers in a loop

Hi everyone, I'm trying to run this code (probably there are better ways to write it, but I'm not a pro at Matlab yet).

params is a vector of 4 parameters, Nstep =365*25, and p_sims is a matrix with Nstep+1 columns and as many raws as the number of simulations (1000 in my case).

I need to integrate a function in two variables, u and s and then use the result to create another matrix.

The problem is that I need to insert in the integrand a vector (the jjj column of p_sims), but this gives me an error (Inner matrix dimensions must agree), whether I use .* to multiply by something containing p_sims, or * How can I solve this?

(I tried with a loop to compute the jjj column of CPrice, but it took ages to run the code, so I was wondering if there was a more efficient way)

for jjj=Nstep+1:-1:1
        CFun = @(s,u) exp(-r*(s-jjj)).*p_sims(:,jjj).^gamma.*floor^(1-gamma)/pi.* ...
            exp(-1i*u.*log(p_sims(:,jjj)/floor)).*exp(1i*params(2).*s.*(-u-1i*gamma));   
CPrice(:,jjj) = X.*real(integral2(CFun,jjj,Nstep+1,0,inf)); 

phi(:,jjj) = p_sims(:,jjj).*(exp((params(2)-r)*(Nstep+1-jjj))*(1-(params(4)+1/2*params(3)^2)/params(1))^ ...
    (-params(1)*(Nstep+1-jjj))-1)  + CPrice(:,jjj);
end

[Walter Roberson]

"inner dimensions must agree" can only occur for * or / and not for .* or ./

[Maria445]

More precisely, with .* the error is "Matrix dimensions must agree", while using * the error is "Inner matrix dimensions must agree"

[Maria445]

This is the screenshot

Answer 1

[Walter Roberson]

The only * operation in your line of code is -r*(s-jjj) . We can see from your code that jjj will be a scalar.

Your s argument is input from the integral2() call. We know from https://www.mathworks.com/help/matlab/ref/integral2.html#inputarg_fun that

"The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations."

so we know that s is a 2D array, of unspecified and variable size.

We are not given any information about the size of r. But clearly since the size of s is variable, if r is not a scalar, there will come a time where size(r,2) ~= size(s,1) . That would lead to the error about inner dimensions not matching. If r were a scalar, there would be no problem at that point, and since there is no other * operation in the line of code, we deduce that whatever r is, it is not a scalar.

[Maria445]

I apologize since I wasn't very clear about the problem.

r is indeed a scalar.

What I meant is that, if I multiply p_sims(:,jjj) using the * operation, I get the error message "Inner matrix dimensions must agree", while, if I use the .* operator (in place of *) to multiply p_sims(:,jjj), what I get is the "Matrix dimensions must agree" error, so I'm stuck either way.

[Maria445]

gamma and floor are scalars as well

[Walter Roberson]

Your s argument is input from the integral2() call. We know from https://www.mathworks.com/help/matlab/ref/integral2.html#inputarg_fun that

"The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations."

so we know that s is a 2D array, of unspecified and variable size.

... So how can you multiply it by a vector p_sims(:,jjj) ?

What is it that you thought you were computing by using the vector p_sims(:,jjj) at that point? If you thought s was a scalar then you would be getting a vector result because of the vector p_sims(:,jjj) and that would mean the result of the function was a vector for a scalar input. Is that your intent? That is, are you trying to integrate a number of different functions at the same location? Because you cannot do that with a single integral2() call.

[Maria445]

You're right, now I see...yes, I was hoping to get a vector as a result, thinking s was a scalar.

Do you have any idea how I could do that? I tried to modify the code like this:

for jjj=Nstep+1:-1:1

   for zzz = 1:size(p_sims,1)
     CFun = @(s,u) exp(-r*(s-jjj))*p_sims(zzz,jjj).^gamma.* ...
       floor^(1-gamma)/pi.*exp(-1i*u*log(p_sims(zzz,jjj)/floor)).* ...
       exp(1i*params(2).*s.*(-u-1i*gamma));
     CPrice(zzz,jjj) = X.*real(integral2(CFun,jjj,Nstep+1,0,inf));
   end
   phi(:,jjj) = p_sims(:,jjj).*(exp((params(2)-r)*(Nstep+1-jjj))* ...
       (1-(params(4)+1/2*params(3)^2)/params(1))^ ...
       (-params(1)*(Nstep+1-jjj))-1)+ CPrice(:,jjj);
end

but it is too computationally expensive, since p_sims is a really big matrix.

[Walter Roberson]

One option that could be experimented with (it might not be any faster at all) would be to convert your 2D integral into a pair of 1D integral() in which the outer integral calls upon the inner one. The inner one could be coded with 'ArrayValued', 'on'. To take care of the fact that the outer integral() will be passed a vector, you might need to arrayfun() the inner integral over its inputs.

... Or possibly you could just calculate your function over a mesh of s and u values, getting out a 3D array (s, u, p_sims column) and then trapz() along the appropriate dimensions. That would not be adaptive, but it might possibly be accurate enough for your purposes.

[Maria445]

Ok, I'm going to try and write the code like you suggested.

Thank you so much for helping me.