From: Vero25 on
Hello,

I have issues to make this code runing. I want to provide a matrix for the variable Rho, so that the result is computed for every element of this matrix.
Can somebody give a hand?

NoRatings=17;
left_bound=-20;
DT(1)=left_bound;
for j=2:NoRatings
DT(j)=norminv(PD(j));
end


NoGrids=17;
rho=0.5; % SHOLD BE A MATRIX!!

NoRows=NoRatings*NoRatings*NoGrids;
result_matrix=zeros(NoRows , 4);
% loop through the ratings and asset-correlation grid points
row_counter=0;
for j=1:NoRatings
for k=1:NoRatings
for m=1 : NoGrids



row_counter=row_counter+1;
% j and k correspnd to RatingA and RatingB (of obligor A and
% obligor B)
RatingA=j;
RatingB=k;


if isequal(j,1) || isequal(k,1)
default_correlation =0;
elseif isequal(rho,1)
PD_joint_default = min(PD(j),PD(k));
default_correlation=(PD_joint_default - PD(j)*PD(k)) / (PD(j)*PD(k)*(1-PD(j))*(1-PD(k)))^0.5;
else
% density is the two-dimensional normal density, given the
% correlation rho
density=['exp(-0.5*(x .* x + y*y -2*x*y*' num2str(rho) ')/(1-' num2str(rho) '*' num2str(rho) ' ))' ] ;

% PD_joint_default is the probability of a joint default (given
% Rating A, Rating B, and rho). The evaluation of the double
% integral with a sufficiently low left boundary in both
% dimensions
PD_joint_default= 1/(2*pi)* 1/(1-rho^2)^0.5 * dblquad(inline(density) ,left_bound,DT(j),left_bound,DT(k));

% resulting correlation of the default event
default_correlation=(PD_joint_default - PD(j)*PD(k)) / (PD(j)*PD(k)*(1-PD(j))*(1-PD(k)))^0.5;
end

result_matrix(row_counter,1) =RatingA;
result_matrix(row_counter,2) =RatingB;
result_matrix(row_counter,3) =rho;
result_matrix(row_counter,4) =default_correlation;

end
end
end