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 | with Generic_Random_Functions,
Generic_Real_Linear_Equations;
with Ada.Exceptions;
with Ada.Strings.Fixed; use Ada.Strings.Fixed;
with Ada.Text_IO;
with Ada.Unchecked_Deallocation;
package body Copulas is
package GRF is new Generic_Random_Functions(Real);
package GRLE is new Generic_Real_Linear_Equations(Real, GRA, Integer_Vector);
use GRF, GRLE;
overriding function Simulate(C: Independent_Copula; gen: Generator) return Real_Vector is
r: Real_Vector(1..C.dim);
begin
-- Build an independent U(0,1) vector
for z in r'Range loop
r(z):= Real(Random(gen));
end loop;
return r;
end Simulate;
overriding function Simulate(C: Independent_Copula; gen: Generator_Vector) return Real_Vector is
r: Real_Vector(1..C.dim);
begin
-- Build an independent U(0,1) vector
for z in r'Range loop
r(z):= Real(Random(gen(z)));
end loop;
return r;
end Simulate;
-- Input : independent U(0,1) vector
-- Output : U(0,1) vector with C's Gaussian dependency.
--
procedure Make_dependent(C: Gauss_Copula; r: in out Real_Vector) is
begin
-- Treat the dimensions having dependencies.
-- 1) Uniform -> Normal
for z in 1..C.dim_dep loop
r(z):= Normal_inverse_CDF(r(z));
end loop;
-- Now r is a draw of an *independant* Normal multivariate random variable.
-- 2) Apply the matrix (Quantitative Risk Management book: Algorithm 3.2)
if C.dim_dep > 0 then
r(1..C.dim_dep):= C.Sqrt_Correl_Matrix.all * r(1..C.dim_dep);
end if;
-- Now r is distributed as a *dependant* Normal
-- 3) Normal -> Uniform
for z in 1..C.dim_dep loop
r(z):= Normal_CDF(r(z));
end loop;
-- Now r is U(0,1) with a Gaussian dependence on 1..dim_dep
end Make_dependent;
overriding function Simulate(C: Gauss_Copula; gen: Generator) return Real_Vector is
r: Real_Vector(1..C.dim);
begin
-- Start with an independent U(0,1) vector
for z in r'Range loop
r(z):= Real(Random(gen));
end loop;
Make_dependent(C, r);
return r;
end Simulate;
overriding function Simulate(C: Gauss_Copula; gen: Generator_Vector) return Real_Vector is
r: Real_Vector(1..C.dim);
begin
for z in r'Range loop
r(z):= Real(Random(gen(z)));
end loop;
Make_dependent(C, r);
return r;
end Simulate;
procedure Dump(A: Real_Matrix; name: String) is
use Ada.Text_IO;
f: File_Type;
begin
Create(f, Out_File, name);
for i in A'Range(1) loop
for j in A'Range(2) loop
Put(f, Real'Image(A(i,j)) & ';');
end loop;
New_Line(f);
end loop;
Close(f);
end Dump;
procedure Construct_as_Gauss(
C : out Copula_access;
dim : Positive;
corr: Real_Matrix
)
is
g: Gauss_Copula(dim);
begin
g.Sqrt_Correl_Matrix:= new Real_Matrix'(Cholesky_Decomposition(corr));
if trace then
Dump(corr, "A.csv");
Dump(g.Sqrt_Correl_Matrix.all, "L.csv");
end if;
g.dim_dep:= corr'Last(1);
--
C:= new Gauss_Copula'(g);
exception
when E : Matrix_Data_Error =>
if Index (Ada.Exceptions.Exception_Message(E), "not positive definite") > 0 then
-- We add an explanation about what it means in terms of correlations.
raise Matrix_Data_Error with
"Matrix is not positive definite: " &
"there is a set of three or more correlations in the matrix that is not consistent. " &
"It is possible to rectify the matrix using the Higham algorithm.";
else
raise; -- Simple re-raise
end if;
end Construct_as_Gauss;
function Get_Cholesky_Matrix(C: Gauss_Copula) return Real_Matrix is
begin
return C.Sqrt_Correl_Matrix.all;
end Get_Cholesky_Matrix;
procedure Dispose_internal is
new Ada.Unchecked_Deallocation(Copula'Class, Copula_access);
procedure Dispose(C: in out Copula_access) is
begin
Dispose_internal(C);
end Dispose;
end Copulas;
|