Code covered by the BSD License

# Credit Risk Modeling with MATLAB

### Ameya Deoras (view profile)

07 Jun 2010 (Updated )

These are the supporting MATLAB files for the MathWorks webinar of the same name.

Credit_VaR.m
%% Credit Risk Analysis
% This script follows a simplified version of the JP Morgan Credit Matrics
% approch.  The MathWorks does not endorse or promote the approach and
% places no guarantees on this code; it is simply provided as an example of
% the type of analyses that are possible.
%
% More complete information about CreditMetrics can be found in many
% places, for instance in the technical document located at:
% http://www.defaultrisk.com/pp_model_20.htm

%% 1. Load data
% The needed data is stored in several places in our spreadsheet.

TransitionMatrix  = dataset( ...
'XLSFile', 'CreditPortfolio.xlsx', ...
'Sheet', 'Transitions and Rates', ...
'Range', 'A1:I8', ...
% This will prove useful:
Ratings = get(TransitionMatrix, 'VarNames');

RecoveryRates = dataset( ...
'XLSFile', 'CreditPortfolio.xlsx', ...
'Sheet', 'Transitions and Rates', ...
'Range', 'A16:C21', ...
% This will prove useful:
Seniorities = get(RecoveryRates, 'ObsNames');

BondData = dataset( ...
'XLSFile', 'CreditPortfolio.xlsx', ...
'Sheet', 'Portfolio Information', ...
% We can optimize our bond data:
BondData.Rating = ordinal(BondData.Rating, [], Ratings);

BondData.Maturity = datenum(BondData.Maturity, 'mm/dd/yyyy');
BondData.Seniority = ordinal(BondData.Seniority, [], Seniorities);

Rates = dataset( ...
'XLSFile', 'CreditPortfolio.xlsx', ...
'Sheet', 'Transitions and Rates', ...
'Range', 'A10:H14', ...

CorrelationMatrix = xlsread( ...
'CreditPortfolio.xlsx', ...
'Correlations', ...
'B2:M13');
% We must fill in the upper half of the correlation matrix.
ctrans = CorrelationMatrix';

CorrelationAlpha = xlsread( ...
'CreditPortfolio.xlsx', ...
'Correlations', ...
'A16');

numScenarios = xlsread( ...
'CreditPortfolio.xlsx', ...
'Risk Analysis', ...
'A7');

[~,valDate] = xlsread( ...
'CreditPortfolio.xlsx', ...
'Risk Analysis', ...
'B7');
valDate = datenum(valDate, 'mm/dd/yyyy');

PercentForVaR = xlsread( ...
'CreditPortfolio.xlsx', ...
'Risk Analysis', ...
'A2:A4');
PercentForVaR = 100 * PercentForVaR;

numAssets = size(BondData,1);

%% 2. Compute original value of portfolio
dates = [valDate; valDate+365; valDate+2*365; valDate+3*365];

% The |prbyzero| function in the Financial Toolbox allows us to value these
% bonds.

for idx = 1 : length(Ratings)-1
ratingsMask = BondData.Rating == Ratings{idx};
BondData.ValuePerBond(ratingsMask) = prbyzero( ...
valDate, Rates.(Ratings{idx}), dates);
end

OriginalPortfolioValue = BondData.ValuePerBond' * BondData.NumberOfBonds;

%% 3. Compute ratings thresholds
% We now need to calculate the ratings thresholds. These represent the
% number of standard deviations an asset price has to move before the
% parent company suffers a rating change.
%
% Using the Asset Value Model, we will compute for each of the ratings the
% Zdefault, Zccc, Zb, Zbb, Zbbb, Za, Zaa, and Zaaa thresholds. We will
%
% $$P_{Default}=P(R<Z_{Default})=\phi(Z_{Default}\mid\sigma)$$
%
% and then iteratively compute the other thresholds using the following
% formula :
%
% $$P_{CCC}=P(Z_{Default}<R<Z_{CCC})=\phi(Z_{CCC}\mid\sigma)-\phi(Z_{Default} \mid\sigma)$$

% Cumulatively add the transition probabilities, starting at the lower end.
ThresholdMatrix = fliplr(cumsum(fliplr(double(TransitionMatrix)), 2));
% Turn these into thresholds.
ThresholdMatrix = norminv(ThresholdMatrix);
% Convert into a dataset array with a similar format as |TransitionMatrix|.
ThresholdMatrix = replacedata(TransitionMatrix, ThresholdMatrix);

%% 4. Generate scenarios
% We will generate a large number of scenarios. To do this, we create a set
% of random numbers correlated by industry within each of the scenarios.
% For each bond we will then add an idiosyncratic random number
% corresponding to the variation within each industry.  The parameter
% governing the relative importance of these two random processes is
% |CorrelationAlpha|, which for the sake of simplicity we have assumed to
% be constant for each bond in our portfolio.

IndustryRandom = mvnrnd( ...
zeros(size(CorrelationMatrix, 1), 1), ...
CorrelationMatrix, ...
numScenarios)';

% By adjusting the weights properly in terms of |CorrelationAlpha|, we
% ensure that the marginal distribution of each company's performance is a
% standard, normally-distributed number.  This mathes the requirements of
% our ratings thresholds.
BondRandom = CorrelationAlpha*IndustryRandom(BondData.Industry,:) + ...
sqrt(1-CorrelationAlpha.^2)*randn(numAssets,numScenarios);

%% 5. Find rating of each bond at the end of the simulations
% We use the threshold levels and the generated scenarios. For each
% scenario, find how the companies fared, and then give them a new rating
% if they have crossed a threshold.
NewRatings = ones(numAssets,numScenarios);

for i = 1 : numAssets
FirmRating = BondData.Rating(i);
Thresholds = [double(ThresholdMatrix(char(FirmRating),:)) -Inf];
for j=2:numel(Thresholds)
RatingIndex = find(BondRandom(i,:) > Thresholds(j) & BondRandom(i,:) < Thresholds(j-1));
NewRatings(i,RatingIndex) = j-1;
end
end

% We no longer need this large variable.
clear BondRandom

NewRatings = ordinal(NewRatings, Ratings);

%% 6. Value bonds based on the simulated ratings
% Now we need to value the bonds for each of the scenarios. This is similar
% to step 2: we price each bond if the bond has not defaulted, and if it
% has, then we calculate how much we expect to recover.
%
% It will be a significant time saver if we pre-calculate the non-default
% prices.  This will save us from having to re-calculate the prices
% multiple times.
NondefaultPrices = prbyzero([BondData.Maturity BondData.Coupon], ...
valDate, double(Rates), dates);

NewPrices = zeros(numAssets,numScenarios);

for Bond = 1 : numAssets
for idx = 1 : length(Ratings)

% Consider nondefault cases:
if idx <= 7;
ratingsMask = NewRatings(Bond,:) == Ratings{idx};

% Consider default cases:
else
Seniority = char(BondData{Bond,'Seniority'});
DistribParameters = double(RecoveryRates(Seniority,{'Mean','STD'}));
defaultMask = NewRatings(Bond,:) == 'Default';
NewPrices(Bond, defaultMask) = ...
100 * betarndms(DistribParameters(1),DistribParameters(2),...
1,NbToDraw);
end
end
end

%% 7. Compute values of simulated portfolios
% The |bsxfun| function is designed to efficiently expand a vector (like
% the 1311x1 "number of bonds" so that it can be multiplied
% element-by-element against a matrix (like the 1311x10000 "New Prices").
NewPrices = bsxfun(@times, NewPrices, BondData.NumberOfBonds);
NewPortfolioValues = sum(NewPrices,1);

%% 8. Display the distribution of possible values for the portfolio

Losses = OriginalPortfolioValue - NewPortfolioValues;

VaR = prctile(Losses, PercentForVaR);
relativeVaR = VaR / OriginalPortfolioValue;

CVaR = zeros(size(VaR));

for i = 1:length(PercentForVaR)
CVaR(i) = mean(Losses(Losses >= VaR(i)));
end

relativeCVaR = CVaR / OriginalPortfolioValue;

figure, hold on,
set(gcf,'Position',[100 100 800 350])
hist(Losses,100);
l1 = plot([VaR(2) VaR(2)], ylim, 'b', 'LineWidth', 2);
l2 = plot([CVaR(2) CVaR(2)], ylim, 'r', 'LineWidth', 2);
set(gca, 'FontSize', 16, 'FontName', 'Calibri', ...
'YTickLabel','')
xlabel('Losses (in USD)')
legend([l1 l2], ...
[int2str(PercentForVaR(2)) '% VaR: ' cur2str(VaR(2),0)], ...
[int2str(PercentForVaR(2)) '% CVaR: ' cur2str(CVaR(2),0)])
hold off