1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
| function results=varin(data)
% data is a column vector of the returns defined as ln(p(t)/p(t-1))*100
% or (p(t)-p(t-1))/p(t-1)*100
% the model is,
% r(t)=m+e(t)
% s(t)^2=v+a*e(t-1)^2+b*s(t-1)^2
% results, is a structure of the constant in mean (m), constant in variance
% (v), arch term (a), garch term (b), the maximum likelihood estimation
% (mle), and a column vector of the conditional standard deviation (sigma).
alpha=[.95;.975;.99;.995;.9975;.999]; % confidence intervals
% matrices for the optimization routine. All the inequalities are
% A*x<=b
Coeff.A=[0 -1 0 0;0 0 -1 0;0 0 0 -1;0 0 1 1;0 0 1 0;0 0 0 1];
Coeff.b=[0;0;0;1;1;1];
% Initialization of the parameter vector
Coeff.c=[0.1;0.1;0.1;0.9];
%mean equation
m=a*(data(t-1)/data(t-1)^2)-d*data(t-1)
% maximum likelihood estimation of Garch(1,1) and Gaussian distribution
function f=garch11nmle(x)
r2=(data-x(1)).^2;
sigma2=zeros(length(data),1);
sigma2(1)=var(data);
for index1=2:length(data)
sigma2(index1)=x(2)+x(3)*r2(index1-1)+x(4)*sigma2(index1-1);
end
f=zeros(length(data),1);
for index2 =1:length(data)
% matlab minimizes the mle so we need an extra minus
f(index2)=-(-1/2*log(2*pi)-1/2.*log(x(2)+x(3)*r2(index2)+x(4)*sigma2(index2))-1/2*r2(index2)/sigma2(index2));
end
f=sum(f);
end
% Optimization
options = optimset('Algorithm','interior-point','Display','off');
[x,fval] = fmincon(@garch11nmle,Coeff.c,Coeff.A,Coeff.b,a,d,[],[],[],[],[],options);
Sigmas=sqrt(sigma2); % standard deviations
% Output structure
results.m=x(1);
results.v=x(2);
results.arch=x(3);
results.garch=x(4);
results.mle=-fval;
results.sigmas=Sigmas;
% Plot of the Value at Risk levels
for index3=1:length(alpha);
varL=x(1)-norminv(alpha(index3)).*Sigmas;
varS=x(1)-norminv(1-alpha(index3)).*Sigmas;
figure
plot(data,'or');
hold on
plot(varL,'-g');
plot(varS,'-b');
hold off
xlabel('Time');
ylabel('VaR and returns');
title(['VaR in sample for alpha=',num2str(alpha(index3),'%6.4f')]);
end
% In sample VaR backtesting
disp([' ','IN SAMPLE VALUE AT RISK BACKTESTING']);
disp('-------------------------------------------------------------------------------------------------------');
disp([' ','SHORT POSITION' ]);
disp('-------------------------------------------------------------------------------------------------------');
disp(['alpha',' ','Success rate',' ','LR value',' ','Kupiec-p',' ','LR Christ',' ','Christ']);
% Short position
for index4=1:length(alpha);
varS=x(1)-norminv(1-alpha(index4)).*Sigmas;
% Kupiec PF proportion of failures test
pval=1-alpha(index4);
Tval=length(data);
xvalS1=(data>varS);
xvalS=sum(xvalS1);
LR_S=-2*(xvalS*log(pval)+(Tval-xvalS)*log(1-pval)-xvalS*log(xvalS/Tval)-(Tval-xvalS)*log(1-xvalS/Tval));
Kupiec_PFS=1-chi2cdf(LR_S,1);
% Christoffersen independence test for conditional coverage
n00=0; n01=0; n10=0; n11=0;
for index5=1:length(xvalS1)-1
if xvalS1(index5)==0 && xvalS1(index5+1)==0
n00=n00+1;
elseif xvalS1(index5)==0 && xvalS1(index5+1)==1
n01=n01+1;
elseif xvalS1(index5)==1 && xvalS1(index5+1)==0
n10=n10+1;
elseif xvalS1(index5)==1 && xvalS1(index5+1)==1;
n11=n11+1;
end
end
p0=n01/(n00+n01);
p1=n11/(n10+n11);
p=(n01+n11)/(n00+n01+n10+n11);
LR_ITS=-2*(n00+n10)*log(1-p)-2*(n01+n11)*log(p)+2*n00*log(1-p0)+2*n01*log(p0)+2*n10*log(1-p1)+2*n11*log(p1);
% LR_ITS=-2*log((1-p)^(n00+n10)*p^(n01+n11)/((1-p0)^n00*p0^n01*(1-p1)^n10*p1^n11));
CIT_S=1-chi2cdf(LR_ITS,1);
fprintf('%6.4f %10.4f %12.4f %14.4e %11.4f %13.4f',alpha(index4),1-xvalS/Tval,LR_S,Kupiec_PFS,LR_ITS,CIT_S);
fprintf('\n');
% disp([num2str(alpha(index4),'%10.4f'),' ',num2str(1-xvalS/Tval,'%10.6f'),' ',num2str(LR_S,'%10.6f'),...
% ' ',num2str(Kupiec_PFS,'%10.6e'),' ',num2str(LR_ITS,'%10.6f'),' ',num2str(CIT_S,'%10.6f')]);
end
disp('-------------------------------------------------------------------------------------------------------');
disp([' ','LONG POSITION' ]);
disp('-------------------------------------------------------------------------------------------------------');
disp(['alpha',' ','Failure rate',' ','LR value',' ','Kupiec-p',' ','LR Christ',' ','Christ']);
% Long position
for index6=1:length(alpha);
varL=x(1)-norminv(alpha(index6)).*Sigmas;
% Kupiec PF proportion of failures test
pval=1-alpha(index6);
Tval=length(data);
xvalL1=(data<varL);
xvalL=sum(xvalL1);
%LR_L=-2*log((pval^xvalL*(1-pval)^(Tval-xvalL))/((xvalL/Tval)^xvalL*(1-xvalL/Tval)^(Tval-xvalL)));
LR_L=-2*(xvalL*log(pval)+(Tval-xvalL)*log(1-pval)-xvalL*log(xvalL/Tval)-(Tval-xvalL)*log(1-xvalL/Tval));
Kupiec_PFL=1-chi2cdf(LR_L,1);
% Christoffersen independence test for conditional coverage
n00=0; n01=0; n10=0; n11=0;
for index7=1:length(xvalL1)-1
if xvalL1(index7)==0 && xvalL1(index7+1)==0
n00=n00+1;
elseif xvalL1(index7)==0 && xvalL1(index7+1)==1
n01=n01+1;
elseif xvalL1(index7)==1 && xvalL1(index7+1)==0
n10=n10+1;
elseif xvalL1(index7)==1 && xvalL1(index7+1)==1;
n11=n11+1;
end
end
p0=n01/(n00+n01);
p1=n11/(n10+n11);
p=(n01+n11)/(n00+n01+n10+n11);
LR_ITL=-2*(n00+n10)*log(1-p)-2*(n01+n11)*log(p)+2*n00*log(1-p0)+2*n01*log(p0)+2*n10*log(1-p1)+2*n11*log(p1);
% LR_ITL=-2*log((1-p)^(n00+n10)*p^(n01+n11)/((1-p0)^n00*p0^n01*(1-p1)^n10*p1^n11));
CIT_L=1-chi2cdf(LR_ITL,1);
fprintf('%6.4f %10.4f %12.4f %14.4e %11.4f %13.4f',alpha(index6),1-xvalL/Tval,LR_L,Kupiec_PFL,LR_ITL,CIT_L);
fprintf('\n');
end
end |
