www.gusucode.com > fininst 案例源码程序 matlab代码 > fininst/Example_CounterpartyCreditRisk.m
%% Counterparty Credit Risk and CVA % % This example shows how to compute the unilateral credit value (valuation) % adjustment (CVA) for a bank holding a portfolio of vanilla interest rate % swaps with several counterparties. CVA is the expected loss on an % over-the-counter contract or portfolio of contracts due to counterparty % default. The CVA for a particular counterparty is defined as the sum % over all points in time of the discounted expected exposure at each % moment multiplied by the probability that the counterparty defaults at % that moment, all multiplied by 1 minus the recovery rate. The CVA % formula is: % % $CVA = (1-R) \int_{0}^{T} discEE(t) dPD(t)$ % % Where |R| is the recovery, |discEE| the discounted expected exposure at % time t, and |PD| the default probability distribution. % % The expected exposure is computed by first simulating many future % scenarios of risk factors for the given contract or portfolio. Risk % factors can be interest rates, as in this example, but will differ based % on the portfolio and can include FX rates, equity or commodity prices, or % anything that will affect the market value of the contracts. Once a % sufficient set of scenarios has been simulated, the contract or portfolio % can be priced on a series of future dates for each scenario. The result % is a matrix, or "cube", of contract values. % % These prices are converted into exposures after taking into account % collateral agreements that the bank might have in place as well as % netting agreements, as in this example, where the values of several % contracts may offset each other, lowering their total exposure. % % The contract values for each scenario are discounted to compute the % discounted exposures. The discounted expected exposures can then be % computed by a simple average of the discounted exposures at each % simulation date. % % Finally, counterparty default probabilities are typically derived from % credit default swap (CDS) market quotes and the CVA for the counterparty % can be computed according to the above formula. We assume that a % counterparty default is independent of its exposure (no wrong-way risk). % % For this example we will work with a portfolio of vanilla interest rate % swaps with the goal of computing the CVA for a particular counterparty. % Copyright 2012-2015 The MathWorks, Inc. %% Read Swap Portfolio % The portfolio of swaps is close to zero value at time t = 0. Each swap % is associated with a counterparty and may or may not be included in a % netting agreement. % Read swaps from spreadsheet swapFile = 'cva-swap-portfolio.xls'; swapData = readtable(swapFile,'Sheet','Swap Portfolio'); swaps = struct( ... 'Counterparty',[], ... 'NettingID',[], ... 'Principal',[], ... 'Maturity',[], ... 'LegRate',[], ... 'LegType',[], ... 'LatestFloatingRate',[], ... 'FloatingResetDates',[]); swaps.Counterparty = swapData.CounterpartyID; swaps.NettingID = swapData.NettingID; swaps.Principal = swapData.Principal; swaps.Maturity = swapData.Maturity; swaps.LegType = [swapData.LegType ~swapData.LegType]; swaps.LegRate = [swapData.LegRateReceiving swapData.LegRatePaying]; swaps.LatestFloatingRate = swapData.LatestFloatingRate; swaps.Period = swapData.Period; swaps.LegReset = ones(size(swaps.LegType)); numSwaps = numel(swaps.Counterparty); %% % For more information on the swap parameters for |Counterparty| and % |NettingID|, see <docid:finance_ug.bu94ri9-1>. For more information on the % swap parameters for |Principal|, |Maturity|, |LegType|, |LegRate|, % |LatestFloatingRate|, |Period|, and |LegReset|, see % <docid:fininst_ug.bu39vqm-1>. %% Create RateSpec from the Interest Rate Curve Settle = datenum('14-Dec-2007'); Tenor = [3 6 12 5*12 7*12 10*12 20*12 30*12]'; ZeroRates = [0.033 0.034 0.035 0.040 0.042 0.044 0.048 0.0475]'; ZeroDates = datemnth(Settle,Tenor); Compounding = 2; Basis = 0; RateSpec = intenvset('StartDates', Settle,'EndDates', ZeroDates, ... 'Rates', ZeroRates,'Compounding',Compounding,'Basis',Basis); figure; plot(ZeroDates, ZeroRates, 'o-'); xlabel('Date'); datetick('keeplimits'); ylabel('Zero rate'); grid on; title('Yield Curve at Settle Date'); %% Set Changeable Simulation Parameters % We can vary here the number of simulated interest rate scenarios we % generate. We set our simulation dates to be more frequent at first, then % turning less frequent further in the future. % Number of Monte Carlo simulations numScenarios = 1000; % Compute monthly simulation dates, then quarterly dates later. simulationDates = datemnth(Settle,0:12); simulationDates = [simulationDates datemnth(simulationDates(end),3:3:74)]'; numDates = numel(simulationDates); %% Compute Floating Reset Dates % For each simulation date, we compute previous floating reset date for % each swap. floatDates = cfdates(Settle-360,swaps.Maturity,swaps.Period); swaps.FloatingResetDates = zeros(numSwaps,numDates); for i = numDates:-1:1 thisDate = simulationDates(i); floatDates(floatDates > thisDate) = 0; swaps.FloatingResetDates(:,i) = max(floatDates,[],2); end %% Setup Hull-White Single Factor Model % The risk factor we will simulate to value our contracts is the zero % curve. For this example we will model the interest rate term structure % using the one-factor Hull-White model. This is a model of the short rate % and is defined as: % % $$ dr = [ \theta(t) - ar ] dt+ \sigma dz $$ % % where % % * $dr$: Change in the short rate after a small change in time, $dt$ % * $a$: Mean reversion rate % * $\sigma$: Volatility of the short rate % * $dz$: A Weiner process (a standard normal process) % * $\theta(t)$: Drift function defined as: % % $$\theta(t) = F_t (0,t) + aF(0,t)+ \frac{\sigma^2}{2a}(1-e^{-2at})$$ % % $F(0,t)$: Instantaneous forward rate at time $t$ % % $F_t (0,t)$: Partial derivative of $F$ with respect to time % % Once we have simulated a path of the short rate we generate a full yield % curve at each simulation date using the formula: % % $$R(t,T) = -\frac{1}{(T-t)} \ln A(t,T) + \frac{1}{(T-t)} B(t,T)r(t)$$ % % $$\ln A(t,T) = \ln \frac{P(0,T)}{P(0,t)} + B(t,T) F(0,t) - \frac{1}{4a^3} \sigma^2 (e^{-aT}-e^{-at} )^2 (e^{2at}-1)$$ % % $$B(t,T) = \frac{1-e^{-a(T-t)}}{a}$$ % % $R(t,T)$: Zero rate at time $t$ for a period of $T-t$ % % $P(t,T)$: Price of a zero coupon bond at time $t$ that pays one dollar at time $T$ % % Each scenario contains the full term structure moving forward through % time, modeled at each of our selected simulation dates. % % Refer to "Calibrating the Hull-White Model Using Market Data" example in % the Financial Instruments Toolbox(TM) Users' Guide for more details on % Hull-White one-factor model calibration. Alpha = 0.2; Sigma = 0.015; hw1 = HullWhite1F(RateSpec,Alpha,Sigma); %% Simulate Scenarios % For each scenario, we simulate the future interest rate curve at each % valuation date using the Hull-White one-factor interest rate model. % Use reproducible random number generator (vary the seed to produce % different random scenarios). prevRNG = rng(0, 'twister'); dt = diff(yearfrac(Settle,simulationDates,1)); nPeriods = numel(dt); scenarios = hw1.simTermStructs(nPeriods, ... 'nTrials',numScenarios, ... 'deltaTime',dt); % Restore random number generator state rng(prevRNG); % Compute the discount factors through each realized interest rate % scenario. dfactors = ones(numDates,numScenarios); for i = 2:numDates tenorDates = datemnth(simulationDates(i-1),Tenor); rateAtNextSimDate = interp1(tenorDates,squeeze(scenarios(i-1,:,:)), ... simulationDates(i),'linear','extrap'); % Compute D(t1,t2) dfactors(i,:) = zero2disc(rateAtNextSimDate, ... repmat(simulationDates(i),1,numScenarios),simulationDates(i-1),-1,3); end dfactors = cumprod(dfactors,1); %% Inspect a Scenario % Create a surface plot of the yield curve evolution for a particular % scenario. i = 20; figure; surf(Tenor, simulationDates, scenarios(:,:,i)) axis tight datetick('y','mmmyy'); xlabel('Tenor (Months)'); ylabel('Observation Date'); zlabel('Rates'); ax = gca; ax.View = [-49 32]; title(sprintf('Scenario %d Yield Curve Evolution\n',i)); %% Compute Mark to Market Swap Prices % For each scenario the swap portfolio is priced at each future simulation % date. Prices are computed using a price approximation function, % |hswapapprox|. It is common in CVA applications to use simplified % approximation functions when pricing contracts due to the performance % requirements of these Monte Carlo simulations. % % Since the simulation dates do not correspond to the swaps cash flow dates % (where the floating rates are reset) we estimate the latest floating rate % with the 1-year rate (all swaps have period 1 year) interpolated between % the nearest simulated rate curves. % % The swap prices are then aggregated into a "cube" of values which % contains all future contract values at each simulation date for each % scenario. The resulting cube of contract prices is a 3 dimensional % matrix where each row represents a simulation date, each column an % contract, and each "page" a different simulated scenario. % Compute all mark-to-market values for this scenario. We use an % approximation function here to improve performance. values = hcomputeMTMValues(swaps,simulationDates,scenarios,Tenor); %% Inspect Scenario Prices % Create a plot of the evolution of all swap prices for a particular % scenario. i = 32; figure; plot(simulationDates, values(:,:,i)); datetick; ylabel('Mark-To-Market Price'); title(sprintf('Swap prices along scenario %d', i)); %% Visualize Simulated Portfolio Values % We plot the total portfolio value for each scenario of our simulation. As % each scenario moves forward in time the values of the contracts will move % up or down depending on how the modeled interest rate term structure % changes. As the swaps get closer to maturity, their values will begin to % approach zero since the aggregate value of all remaining cash flows will % decrease after each cash flow date. % View portfolio value over time figure; totalPortValues = squeeze(sum(values, 2)); plot(simulationDates,totalPortValues); title('Total MTM Portfolio Value for All Scenarios'); datetick('x','mmmyy') ylabel('Portfolio Value ($)') xlabel('Simulation Dates') %% Compute Exposure by Counterparty % The exposure of a particular contract (i) at time t is the maximum of the % contract value (Vi) and 0: % % $$ E_i (t) = \max \{ V_i (t),0 \} $$ % % And the exposure for a particular counterparty is simply a sum of the % individual contract exposures: % % $$ E_{cp}(t) = \sum E_i (t) = \sum \max \{ V_i (t),0 \} $$ % % In the presence of netting agreements, however, contracts are aggregated % together and can offset each other. Therefore the total exposure of all % contracts in a netting agreement is: % % $$ E_{na}(t) = \max \{ \sum V_i (t), 0 \} $$ % % We compute these exposures for the entire portfolio as well as each % counterparty at each simulation date using the |creditexposures| % function. % % Unnetted contracts are indicated using a NaN in the NettingID vector. % Exposure of an unnetted contract is equal to the market value of the % contract if it has positive value, otherwise it is zero. % % Contracts included in a netting agreement have their values aggregated % together and can offset each other. See the references for more details % on computing exposure from mark-to-market contract values. [exposures, expcpty] = creditexposures(values,swaps.Counterparty, ... 'NettingID',swaps.NettingID); %% % We plot the total portfolio exposure for each scenario in our simulation. % Similar to the plot of contract values, the exposures for each scenario % will approach zero as the swaps mature. % View portfolio exposure over time figure; totalPortExposure = squeeze(sum(exposures,2)); plot(simulationDates,totalPortExposure); title('Portfolio Exposure for All Scenarios'); datetick('x','mmmyy') ylabel('Exposure ($)') xlabel('Simulation Dates') %% Exposure Profiles % Several exposure profiles are useful when analyzing the potential future % exposure of a bank to a counterparty. Here we compute several % (non-discounted) exposure profiles per counterparty as well as for the % entire portfolio. % % * |PFE| : Potential Future Exposure : A high percentile (95%) of the % distribution of exposures at any particular future date. Also called % Peak Exposure (PE) % * |MPFE| : Maximum Potential Future Exposure : The maximum PFE across all % dates % * |EE| : Expected Exposure : The mean (average) of the distribution of % exposures at each date % * |EPE| : Expected Positive Exposure : Weighted average over time of the % expected exposure % * |EffEE| : Effective Expected Exposure : The maximum expected exposure at % any time, t, or previous time % * |EffEPE| : Effective Expected Positive Exposure : The weighted average % of the effective expected exposure % % For further definitions, see for example Basel II document in references. % Compute entire portfolio exposure portExposures = sum(exposures,2); % Compute exposure profiles for each counterparty and entire portfolio cpProfiles = exposureprofiles(simulationDates,exposures); portProfiles = exposureprofiles(simulationDates,portExposures); %% % We visualize the exposure profiles, first for the entire portfolio, then % for a particular counterparty. % Visualize portfolio exposure profiles figure; plot(simulationDates,portProfiles.PFE, ... simulationDates,portProfiles.MPFE * ones(numDates,1), ... simulationDates,portProfiles.EE, ... simulationDates,portProfiles.EPE * ones(numDates,1), ... simulationDates,portProfiles.EffEE, ... simulationDates,portProfiles.EffEPE * ones(numDates,1)); legend({'PFE (95%)','Max PFE','Exp Exposure (EE)','Time-Avg EE (EPE)', ... 'Max past EE (EffEE)','Time-Avg EffEE (EffEPE)'}) datetick('x','mmmyy') title('Portfolio Exposure Profiles'); ylabel('Exposure ($)') xlabel('Simulation Dates') %% % Visualize exposure profiles for a particular counterparty cpIdx = find(expcpty == 5); figure; plot(simulationDates,cpProfiles(cpIdx).PFE, ... simulationDates,cpProfiles(cpIdx).MPFE * ones(numDates,1), ... simulationDates,cpProfiles(cpIdx).EE, ... simulationDates,cpProfiles(cpIdx).EPE * ones(numDates,1), ... simulationDates,cpProfiles(cpIdx).EffEE, ... simulationDates,cpProfiles(cpIdx).EffEPE * ones(numDates,1)); legend({'PFE (95%)','Max PFE','Exp Exposure (EE)','Time-Avg EE (EPE)', ... 'Max past EE (EffEE)','Time-Avg EffEE (EffEPE)'}) datetick('x','mmmyy','keeplimits') title(sprintf('Counterparty %d Exposure Profiles',cpIdx)); ylabel('Exposure ($)') xlabel('Simulation Dates') %% Discounted Exposures % We compute the discounted expected exposures using the discount factors % from each simulated interest rate scenario. The discount factor for a % given valuation date in a given scenario is the product of the % incremental discount factors from one simulation date to the next, along % the interest rate path of that scenario. % Get discounted exposures per counterparty, for each scenario discExp = zeros(size(exposures)); for i = 1:numScenarios discExp(:,:,i) = bsxfun(@times,dfactors(:,i),exposures(:,:,i)); end % Discounted expected exposure discProfiles = exposureprofiles(simulationDates,discExp, ... 'ProfileSpec','EE'); %% % We plot the discounted expected exposures for the aggregate portfolio as % well as for each counterparty. % Aggregate the discounted EE for each counterparty into a matrix discEE = [discProfiles.EE]; % Portfolio discounted EE figure; plot(simulationDates,sum(discEE,2)) datetick('x','mmmyy','keeplimits') title('Discounted Expected Exposure for Portfolio'); ylabel('Discounted Exposure ($)') xlabel('Simulation Dates') % Counterparty discounted EE figure; plot(simulationDates,discEE) datetick('x','mmmyy','keeplimits') title('Discounted Expected Exposure for Each Counterparty'); ylabel('Discounted Exposure ($)') xlabel('Simulation Dates') %% Calibrating Probability of Default Curve for Each Counterparty % The default probability for a given counterparty is implied by the % current market spreads of the counterparty's CDS. We use the function % |cdsbootstrap| to generate the cumulative probability of default at each % simulation date. % Import CDS market information for each counterparty CDS = readtable(swapFile,'Sheet','CDS Spreads'); disp(CDS); CDSDates = datenum(CDS.Date); CDSSpreads = table2array(CDS(:,2:end)); ZeroData = [RateSpec.EndDates RateSpec.Rates]; % Calibrate default probabilities for each counterparty DefProb = zeros(length(simulationDates), size(CDSSpreads,2)); for i = 1:size(DefProb,2) probData = cdsbootstrap(ZeroData, [CDSDates CDSSpreads(:,i)], ... Settle, 'probDates', simulationDates); DefProb(:,i) = probData(:,2); end % We plot of the cumulative probability of default for each counterparty. figure; plot(simulationDates,DefProb) title('Default Probability Curve for Each Counterparty'); xlabel('Date'); grid on; ylabel('Cumulative Probability') datetick('x','mmmyy') ylabel('Probability of Default') xlabel('Simulation Dates') %% CVA Computation % % The Credit Value (Valuation) Adjustment (CVA) formula is: % % $CVA = (1-R) \int_{0}^{T} discEE(t) dPD(t)$ % % Where |R| is the recovery, |discEE| the discounted expected exposure at % time t, and |PD| the default probability distribution. This assumes the % exposure is independent of default (no wrong-way risk), and it also % assumes the exposures were obtained using risk-neutral probabilities. % % Here we approximate the integral with a finite sum over the valuation % dates as: % % $CVA (approx) = (1-R) \sum_{i=2}^n discEE(t_i) (PD(t_i)-PD(t_{i-1}))$ % % where t_1 is todays date, t_2, ...,t_n the future valuation dates. % % We assume CDS info corresponds to counterparty with index cpIdx. The % computed CVA is the present market value of our credit exposure to % counterparty cpIdx. For this example we set the recovery rate at 40%. Recovery = 0.4; CVA = (1-Recovery) * sum(discEE(2:end,:) .* diff(DefProb)); for i = 1:numel(CVA) fprintf('CVA for counterparty %d = $%.2f\n',i,CVA(i)); end figure; bar(CVA); title('CVA for each counterparty'); xlabel('Counterparty'); ylabel('CVA $'); grid on; %% References % % # Pykhtin, Michael, and Steven Zhu, _A Guide to Modeling Counterparty % Credit Risk_, GARP, July/August 2007, issue 37, pp. 16-22. % # Pykhtin, Michael, and Dan Rosen, _Pricing Counterparty Risk at the % Trade Level and CVA_, 2010. % # Basel II: http://www.bis.org/publ/bcbs128.pdf page 256