Unable to perform numerical integration when matrices are involved.

I have been badly stuck up with the following code.

if true
  n=35; 
T=15; 
g=.111;
K=1/g;
d=2; 
rbar=0.04;
mubar=0; 
R=[1, 0; 0, 0];
M=[0, 0; 0, 1];
X012=0.002;
X0=[0.01, X012; X012, 0.001];
H=[-0.5, 0.4; 0.007, -0.008];
Q=[0.06 -0.0006; -0.06, 0.006];
Scheck7=0;
beta=3;
SIGMA=[0.006811791307233, -0.000407806990090; -0.000407806990090, 0.00039291243623];
PSICURL=[0.011526035149236, 0.758273970303934; 0.013935191202735, 0.955423605940771];   
%Specifications of the arrays for the Levy Assumption
S0i=zeros(n,1);
S0i2=zeros(n,1);
Yn0=0; %Used in the first lower bound
Yi=zeros(n,1);
%Arrays for the MC estimate
S1MC=zeros(n,1);
POW=zeros(n,1);
SPOW=0;
SPHI=0;
SPSI=0;
% Parameters used for the new approach
% To define the Matrices PSIi in an array
delta=0.75; %The damping factor
PSIi = cell(1, n);
PHIi=zeros(n,1);
% Prameters for the symbolic sum
% Remember our x is Gamma1
%syms x ik
% Computation of Rho
Num=(Q(1,1)*Q(1,2)+Q(2,2)*Q(2,1))*X0(1,2);
Denom=sqrt((Q(1,1)^2+Q(2,1)^2)*X0(1,1)*(Q(2,2)^2+Q(1,2)^2)*X0(2,2));
Rho=Num/Denom
% Specification of new MATRICES Aij's
ANEW=expm([T.*H,T.*(2*(Q'*Q));T.*(R+M),T.*(-(H'))]);
A11=ANEW(1:2,1:2); %That is good
A21=ANEW(3:4,1:2);
A12=ANEW(1:2,3:4);%intersection
A22=ANEW(3:4,3:4);
% Computationof psi(T) and phi(T) (Otherwise with varying i;run in a loop)
C22=inv(A22);
PSIT=C22*A21;
PHIT=beta*(log(det(A22))+T*trace(H'))/2;
% Computation of SZCB's Pcurl(0,T)
C=trace(PSIT*X0);
SZCB=exp(-(rbar+mubar)*T)*exp(-PHIT-C)
SIGMAi=inv(SIGMA);
THETA1=SIGMAi*(PSICURL'*X0*PSICURL);
for i=2:n;
       APHNEW=expm([(i-1).*H,(i-1).*(2*(Q'*Q));(i-1).*(R+M),(i-1).*(-(H'))]);
       APHNEW11=APHNEW(1:2,1:2); %That is good
       APHNEW21=APHNEW(3:4,1:2);
       APHNEW12=APHNEW(1:2,3:4);
       APHNEW22=APHNEW(3:4,3:4);
       % Computation of PSIi and PHIi  
       BPHNEW22=inv(APHNEW22);
       PSIi{i}=APHNEW22\APHNEW21;
       M7=PSIi{i};
       SPSI=SPSI+PSIi{i};
       PHIi(i)=beta*(log(det(APHNEW22))+(i-1)*trace(H'))/2; 
       S0i(i)=exp(-((rbar+mubar)*(i-1)+PHIi(i)));
       S0i2(i)=(-((rbar+mubar)*(i-1)+PHIi(i))); 
       Yn0=Yn0+S0i2(i);
end
SIGMAi=inv(SIGMA);
THETA1=SIGMAi*(PSICURL'*X0*PSICURL);
fun4 = @(Gamma1) exp(-(1i.*Gamma1).*log(K-1)).*exp((1i.*Gamma1+(delta+1)).*Yn0).*(exp(trace(1i.*(THETA1*inv(eye(2)-2i.*SIGMA*((-(Gamma1-1i.*(delta+1))).*SPSI))*SIGMA*((-(Gamma1-1i.*(delta+1))).*SPSI))))./((det(eye(2)-2i.*SIGMA*((-(Gamma1-1i.*(delta+1))).*SPSI))).^(beta./2)))./(delta.^2+delta-Gamma1.^2+(1i.*Gamma1.*(2.*delta+1)));
qty = quadgk(fun4,0,inf)
qty1= exp(-(delta*log(K-1)))*qty/pi
LB1=g*SZCB*qty1
format 'long'
end

The numerical integration involves matrices. I am unable to use quadgk. Can anyone please help?

Thanks in advance.

Answer 1

[Torsten]

Try whether this code works for you:

function main
qty = quadgk(@fun4,0,inf)

function y=fun4(gamma1vec)
<<Put everyting of your code here up to the line THETA1=SIGMAi*(PSICURL'*X0*PSICURL);>> 
for i=1:numel(gamma1vec)
  gamma1=gamma1vec(i);
  y(i)=exp(-(1i.*gamma1).*log(K-1)).*exp((1i.*gamma1+(delta+1)).*Yn0).*(exp(trace(1i.*(THETA1*inv(eye(2)-2i.*SIGMA*((-(gamma1-1i.*(delta+1))).*SPSI))*SIGMA*((-(gamma1-1i.*(delta+1))).*SPSI))))./((det(eye(2)-2i.*SIGMA*((-(gamma1-1i.*(delta+1))).*SPSI))).^(beta./2)))./(delta.^2+delta-gamma1.^2+(1i.*gamma1.*(2.*delta+1)));
end

Best wishes

Torsten.

[RB]

Thank you so much. It does work. I completely accept your answer. I appreciate your efforts. Kind regards, RB

[Torsten]

function main
qty = quadgk(@fun4,0,inf)
[y, delta, K, g, SZCB] = fun4(1);
qty1 = exp(-(delta*log(K-1)))*qty/pi;
LB1=g*SZCB*qty1;

function [y, delta, K, g, SZCB] = fun4(gamma1vec)
<< Include everything up to the end of the for-loop >>

Best wishes

Torsten.

[RB]

Thank you once again.

Regards, RB