www.gusucode.com > risk 案例源码程序 matlab代码 > risk/ValueatRiskEstimationandBacktestingExample.m

    %% Value-at-Risk Estimation and Backtesting 
%
% This example shows how to estimate the value-at-risk (VaR) using three 
% methods, and how to perform a VaR backtesting analysis. The three methods
% are:
%
% # Normal distribution 
% # Historical simulation
% # Exponential weighted moving average (EWMA) 
% 
% Value-at-risk is a statistical method that quantifies the risk level
% associated with a portfolio. The VaR measures the maximum amount of 
% loss over a specified time horizon and at a given confidence level. 
%
% Backtesting measures the accuracy the VaR calculations. Using VaR methods, the loss 
% forecast is calculated and then compared to the actual losses at the 
% end of the next day. The degree of difference between the predicted and actual losses 
% indicates whether the VaR model is underestimating or overestimating the 
% risk. As such, backtesting looks retrospectively at data and helps to assess 
% the VaR model.
% 
% The three estimation methods used in this example estimate the VaR at 95% and 99% 
% confidence levels.

%% Load the Data and Define the Test Window
% Load the data. The data used in this example is from a time series of returns on the S&P index
% between 1993 and 2003.

%%
load VaRExampleData.mat
Returns = diff(sp)./sp(1:end-1);
DateReturns = dates(2:end);
SampleSize = length(Returns); 

%%
% Define the estimation window as 250 trading days. The test window starts 
% on the first day in 1996 and runs through the end of the sample.

%%
TestWindowStart      = find(year(DateReturns)==1996,1);
TestWindow           = TestWindowStart : SampleSize;
EstimationWindowSize = 250; 

%%
% For a VaR confidence level of 95% and 99%, set the complement of the VaR level.
pVaR = [0.05 0.01];

%%
% These values mean that there is at most a 5% and 1% probability, respectively,
% that the loss incurred will be greater than the maximum threshold 
% (that is, greater than the VaR).

%% Compute the VaR Using the Normal Distribution Method
% For the normal distribution method, assume that the profit and loss of the portfolio
% is normally distributed. Using this assumption, compute the VaR by 
% multiplying the _z_-score, at each confidence level by the standard
% deviation of the returns.
% Because VaR backtesting looks retrospectively at data, the VaR
% "today" is computed based on values of the returns in the last _N_ = 250
% days leading to, but not including, "today."

Zscore   = norminv(pVaR);
Normal95 = zeros(length(TestWindow),1);
Normal99 = zeros(length(TestWindow),1);

for t = TestWindow
    i = t - TestWindowStart + 1;
    EstimationWindow = t-EstimationWindowSize:t-1;
    Sigma = std(Returns(EstimationWindow)); 
    Normal95(i) = -Zscore(1)*Sigma;
    Normal99(i) = -Zscore(2)*Sigma;
end

figure;
plot(DateReturns(TestWindow),[Normal95 Normal99])
xlabel('Date')
ylabel('VaR')
legend({'95% Confidence Level','99% Confidence Level'},'Location','Best')
title('VaR Estimation Using the Normal Distribution Method')

%%
% The normal distribution method is also known as parametric VaR because 
% its estimation involves computing a parameter for the standard 
% deviation of the returns. The advantage of the normal distribution method
% is its simplicity. However, the weakness of the normal distribution
% method is the assumption that returns are
% normally distributed. Another name for the normal distribution method is the variance-covariance approach.

%% Compute the VaR Using the Historical Simulation Method
% Unlike the normal distribution method, the historical simulation (HS) is
% a nonparametric method. It does not assume a particular 
% distribution of the asset returns. Historical simulation forecasts risk
% by assuming that past profits and losses can be used as the distribution
% of profits and losses for the next period of returns. The VaR "today" is computed as the 
% _p_ th-quantile of the last  _N_  returns prior to "today."

Historical95 = zeros(length(TestWindow),1);
Historical99 = zeros(length(TestWindow),1);

for t = TestWindow
    i = t - TestWindowStart + 1;
    EstimationWindow = t-EstimationWindowSize:t-1;
    X = Returns(EstimationWindow);
    Historical95(i) = -quantile(X,pVaR(1)); 
    Historical99(i) = -quantile(X,pVaR(2)); 
end

figure;
plot(DateReturns(TestWindow),[Historical95 Historical99])
ylabel('VaR')
xlabel('Date')
legend({'95% Confidence Level','99% Confidence Level'},'Location','Best')
title('VaR Estimation Using the Historical Simulation Method')

%%
% The preceding figure shows that the historical simulation curve has a 
% piecewise constant profile. The reason for this is that quantiles do not change 
% for several days until extreme events occur. Thus, the historical simulation method is slow to react
% to changes in volatility.

%% Compute the VaR Using the Exponential Weighted Moving Average Method (EWMA)
% The first two VaR methods assume that all past returns carry the 
% same weight. The exponential weighted moving average (EWMA) method assigns
% nonequal weights, particularly exponentially decreasing weights. The 
% most recent returns have higher weights because they influence "today's" 
% return more heavily than returns further in the past.
% The formula for the EWMA variance over an estimation window of size $W_E$
% is:
%
% $$\hat{\sigma}^2_t=\frac{1}{c}\sum_{i=1}^{W_E}\lambda^{i-1} y^2_{t-i}$$
%
% where $c$ is a normalizing constant:
%
% $$c=\sum_{i=1}^{W_E}\lambda^{i-1} = \frac{1-\lambda^{W_E}}{1-\lambda}\quad\rightarrow \frac{1}{1-\lambda} ~ as ~ W_E\rightarrow\infty$$
%
% For convenience, we assume an infinitely large estimation window to
% approximate the variance:
%
% $$\hat{\sigma}^2_t\approx(1-\lambda)(y^2_{t-1}+\sum^{\infty}_{i=2}\lambda^{i-1}y^2_{t-i})=(1-\lambda)y^2_{t-1}+\lambda\hat{\sigma}^2_{t-1}$$
%
% A value of the decay factor frequently used in practice is 0.94. This is 
% the value used in this example. For more information, see References.
%
%%
% Initiate the EWMA using a warm-up phase to set up the standard
% deviation.
Lambda = 0.94;
Sigma2     = zeros(length(Returns),1);
Sigma2(1)  = Returns(1)^2;

for i = 2 : (TestWindowStart-1)
    Sigma2(i) = (1-Lambda) * Returns(i-1)^2 + Lambda * Sigma2(i-1);
end

%% 
% Use the EWMA in the test window to estimate the VaR.
Zscore = norminv(pVaR);
EWMA95 = zeros(length(TestWindow),1);
EWMA99 = zeros(length(TestWindow),1);

for t = TestWindow 
    k     = t - TestWindowStart + 1;
    Sigma2(t) = (1-Lambda) * Returns(t-1)^2 + Lambda * Sigma2(t-1);
    Sigma = sqrt(Sigma2(t));
    EWMA95(k) = -Zscore(1)*Sigma;
    EWMA99(k) = -Zscore(2)*Sigma;
end

figure;
plot(DateReturns(TestWindow),[EWMA95 EWMA99])
ylabel('VaR')
xlabel('Date')
legend({'95% Confidence Level','99% Confidence Level'},'Location','Best')
title('VaR Estimation Using the EWMA Method')

%%
% In the preceding figure, the EWMA reacts very quickly to
% periods of large (or small) returns.
%
%% VaR Backtesting
% In the first part of this example, VaR was estimated over the test window
% with three different methods and at two different VaR confidence levels. 
% The goal of VaR backtesting is to evaluate the performance of VaR models. 
% A VaR estimate at 95% confidence is violated only about 5% of the 
% time, and VaR failures do not cluster. Clustering of VaR failures 
% indicates the lack of independence across time because the VaR models are slow to react 
% to changing market conditions.
%
% A common first step in VaR backtesting analysis is to plot the returns and
% the VaR estimates together. Plot all three methods at the 95% confidence level and
% compare them to the returns.

%%
ReturnsTest = Returns(TestWindow);
DatesTest   = DateReturns(TestWindow);
figure;
plot(DatesTest,[ReturnsTest -Normal95 -Historical95 -EWMA95])
ylabel('VaR')
xlabel('Date')
legend({'Returns','Normal','Historical','EWMA'},'Location','Best')
title('Comparison of returns and VaR at 95% for different models')

%%
% To highlight how the different approaches react differently to changing 
% market conditions, you can zoom in on the time series where there is a 
% large and sudden change in the value of returns. For example, around 
% August 1998:

ZoomInd   = (DatesTest >= '5-Aug-1998') & (DatesTest <= '31-Oct-1998');
VaRData   = [-Normal95(ZoomInd) -Historical95(ZoomInd) -EWMA95(ZoomInd)];
VaRFormat = {'-','--','-.'};
figure;
bar(datenum(DatesTest(ZoomInd)),ReturnsTest(ZoomInd),'FaceColor',[0.6 0.6 0.6]);
hold on
for i = 1 : size(VaRData,2)
    stairs(datenum(DatesTest(ZoomInd))-0.5,VaRData(:,i),VaRFormat{i});
end
ylabel('VaR')
xlabel('Date')
legend({'Returns','Normal','Historical','EWMA'},'Location','Best','AutoUpdate','Off')
title('95% VaR violations for different models')
datetick('x','keeplimits');

%%
% A VaR failure or violation happens when the returns have a negative VaR.
% A closer look around August 27th to August 31st shows a
% significant dip in the returns. On the dates starting from August 27th 
% onward, the EWMA follows the trend of the returns closely and 
% more accurately.  Consequently, EWMA has fewer VaR violations (two) 
% compared to the normal distribution approach (seven violations) or the 
% historical simulation method (eight violations). 
%
%% 
% Besides visual tools, you can use statistical tests for VaR backtesting.
% In Risk Management Toolbox(TM), a |varbacktest| object supports
% multiple statistical tests for VaR backtesting analysis.
% In this example, start by comparing the different test results for the 
% normal distribution approach at the 95% and 99% VaR levels.
vbt = varbacktest(ReturnsTest,[Normal95 Normal99],'PortfolioID','S&P','VaRID',...
    {'Normal95','Normal99'},'VaRLevel',[0.95 0.99]);
summary(vbt)
%%
% The summary report shows that the observed level is close enough to the defined
% VaR level. The 95% and 99% VaR levels have at most 
% |(1-VaR_level) x _N_|   expected failures, where _N_ is the number of 
% observations. The failure ratio shows that the |Normal95| VaR level is within range, 
% whereas the |Normal99| VaR Level is imprecise and under-forecasts the risk.
% To run all tests supported in |varbacktest|, use |runtests|.
%
runtests(vbt)

%%
% The 95% VaR passes the frequency tests, such as traffic light, binomial 
% and proportion of failures tests (|TL|, |Bin|, and |POF| columns. The 99% 
% VaR does not pass these same tests, as indicated by the |yellow| and
% |reject| results. Both confidence levels got rejected in the conditional
% coverage independence, and time between failures independence (|CCI| and
% |TBFI| columns. This result suggests that the VaR violations are not
% independent, and there are probably periods with multiple failures in a
% short span. Also, one failure may make it more likely that other failures will follow in subsequent 
% days. For more information on the tests methodologies and the interpretation 
% of results, see <docid:risk_ug.bvaa3t4> and the individual tests.
%
% Using a |varbacktest| object, run the same tests on the portfolio for
% the three approaches at both VaR confidence levels.

vbt = varbacktest(ReturnsTest,[Normal95 Historical95 EWMA95 Normal99 Historical99 ...
    EWMA99],'PortfolioID','S&P','VaRID',{'Normal95','Historical95','EWMA95',...
    'Normal99','Historical99','EWMA99'},'VaRLevel',[0.95 0.95 0.95 0.99 0.99 0.99]);
runtests(vbt)

%%
% The results are similar to the previous results, and at the 95% level, the frequency
% results are generally acceptable. However, the frequency results at the 
% 99% level are generally rejections. Regarding independence, most tests 
% pass the conditional coverage independence test (<docid:risk_ug.bvabiyt-1>), which tests for 
% independence on consecutive days. Notice that all tests fail the time between 
% failures independence test (<docid:risk_ug.bvabi29-1>), which takes into account the times 
% between all failures. This result suggests that all methods have issues with the 
% independence assumption.
%
% To better understand how these results change given market conditions, 
% look at the years 2000 and 2002 for the 95% VaR confidence level.
%
Ind2000 = (year(DatesTest) == 2000);
vbt2000 = varbacktest(ReturnsTest(Ind2000),[Normal95(Ind2000) Historical95(Ind2000) EWMA95(Ind2000)],...
   'PortfolioID','S&P, 2000','VaRID',{'Normal','Historical','EWMA'});
runtests(vbt2000)

%%
Ind2002 = (year(DatesTest) == 2002);
vbt2002 = varbacktest(ReturnsTest(Ind2002),[Normal95(Ind2002) Historical95(Ind2002) EWMA95(Ind2002)],...
   'PortfolioID','S&P, 2002','VaRID',{'Normal','Historical','EWMA'});
runtests(vbt2002)

%%
% For the year 2000, all three methods pass all the tests. However, for the 
% year 2002, the test results are mostly rejections for all methods. The EWMA 
% method seems to perform better in 2002, yet all methods fail the 
% independence tests.
%
% To get more insight into the independence tests, look into the
% conditional coverage independence (<docid:risk_ug.bvabiyt-1>) and the time between failures 
% independence (<docid:risk_ug.bvabi29-1>) test details for the year 2002.
% To access the test details for all tests, run the individual
% test functions.
%%
cci(vbt2002)
%%
% In the CCI test, the probability _p_ |01| of having a failure at 
% time _t_, knowing that there was no failure at time _t_-1 is given by
%
% $$ p_{01} = \frac{N_{01}}{N_{01}+N_{00}}$$
%
% The probability _p_ |11| of having a failure at time _t_, knowing 
% that there was failure at time _t_-1 is given by
%
% $$ p_{11} = \frac{N_{11}}{N_{11}+N_{10}}$$
%
% From the |N00|, |N10|, |N01|, |N11| columns in the test results, the value of 
% _p_ |01| is at around 5% for the three methods, yet the values of _p_ |11| are 
% above 20%. Because there is evidence that a failure is followed by another 
% failure much more frequently than 5% of the time, this 
% CCI test fails.
%
% In the time between failures independence test, look at the minimum, 
% maximum, andquartiles of the distribution of times between failures, 
% in the columns |TBFMin|, |TBFQ1|, |TBFQ2|, |TBFQ3|, |TBFMax|. 
%
%%
tbfi(vbt2002)

%%
% For a VaR level of 95%, you expect an average time between failures 
% of 20 days, or one failure every 20 days. However, the median of the time 
% between failures for the year 2002 ranges between 5 and 7.5 for the three 
% methods. This result suggests that half of the time, two consecutive failures occur within 
% 5 to 7 days, much more frequently than the 20 expected days. Consequently, 
% more test failures occur. For the normal method, the first quartile is 1, 
% meaning that 25% of the failures occur on consecutive days.
%
%% References
% Nieppola, O. _Backtesting Value-at-Risk Models_. Helsinki School of
% Economics. 2009.
%
% Danielsson, J. _Financial Risk Forecasting: The Theory and Practive of Forecasting Market 
% Risk, with Implementation in R and MATLAB(R)_. Wiley 
% Finance, 2012.