1 | /*---------------------------------------------------------------------------- |
---|
2 | ADOL-C -- Automatic Differentiation by Overloading in C++ |
---|
3 | File: liborpar.cpp |
---|
4 | Revision: $Id: liborpar.cpp 708 2016-07-12 08:18:44Z kulshres $ |
---|
5 | Contents: example for differentiation of OpemMP parallel programs |
---|
6 | parallel version |
---|
7 | |
---|
8 | Copyright (c) Andrea Walther |
---|
9 | |
---|
10 | This file is part of ADOL-C. This software is provided as open source. |
---|
11 | Any use, reproduction, or distribution of the software constitutes |
---|
12 | recipient's acceptance of the terms of the accompanying license file. |
---|
13 | |
---|
14 | ---------------------------------------------------------------------------*/ |
---|
15 | |
---|
16 | /* Program to compute deltas and vegas of swaption portfolio |
---|
17 | from forward and reverse mode pathwise sensitivities |
---|
18 | in parallel written by Andrea Walther in 2008-11 based on |
---|
19 | code written by written by Mike Giles in 2005-7 which is |
---|
20 | again based on code written by Zhao and Glasserman at |
---|
21 | Columbia University */ |
---|
22 | |
---|
23 | using namespace std; |
---|
24 | |
---|
25 | #include <stdlib.h> |
---|
26 | #include <stdio.h> |
---|
27 | #include <ctime> |
---|
28 | #include <cmath> |
---|
29 | |
---|
30 | #include "adolc/adolc.h" |
---|
31 | |
---|
32 | #ifdef _OPENMP |
---|
33 | #include <omp.h> |
---|
34 | #include <adolc/adolc_openmp.h> |
---|
35 | #endif |
---|
36 | |
---|
37 | /* calculate path values */ |
---|
38 | |
---|
39 | template <typename ADdouble> |
---|
40 | void path_calc(const int N, const int Nmat, const double delta, |
---|
41 | ADdouble L[], const double lambda[], ADdouble z[]) |
---|
42 | { |
---|
43 | int i, n; |
---|
44 | double lam, con1; |
---|
45 | ADdouble v, vrat; |
---|
46 | ADdouble sqez; |
---|
47 | |
---|
48 | for(n=0; n<Nmat; n++) { |
---|
49 | sqez = sqrt(delta)*z[n]; |
---|
50 | |
---|
51 | v = 0.0; |
---|
52 | for (i=n+1; i<N; i++) { |
---|
53 | lam = lambda[i-n-1]; |
---|
54 | con1 = delta*lam; |
---|
55 | v += (con1*L[i])/(1.0+delta*L[i]); |
---|
56 | vrat = exp(con1*v + lam*(sqez-0.5*con1)); |
---|
57 | L[i] = L[i]*vrat; |
---|
58 | } |
---|
59 | } |
---|
60 | } |
---|
61 | |
---|
62 | /* calculate the portfolio value v */ |
---|
63 | |
---|
64 | template <typename ADdouble> |
---|
65 | void portfolio(const int N, const int Nmat, const double delta, |
---|
66 | const int Nopt, const int maturities[], |
---|
67 | const double swaprates[], |
---|
68 | const ADdouble L[], ADdouble& v ) |
---|
69 | { |
---|
70 | int i, m, n; |
---|
71 | ADdouble b, s, swapval, *B, *S; |
---|
72 | B = new ADdouble[N]; |
---|
73 | S = new ADdouble[N]; |
---|
74 | |
---|
75 | b = 1.0; |
---|
76 | s = 0.0; |
---|
77 | |
---|
78 | for (n=Nmat; n<N; n++) { |
---|
79 | b = b/(1.0+delta*L[n]); |
---|
80 | s = s + delta*b; |
---|
81 | B[n] = b; |
---|
82 | S[n] = s; |
---|
83 | } |
---|
84 | |
---|
85 | v = 0; |
---|
86 | |
---|
87 | for (i=0; i<Nopt; i++){ |
---|
88 | m = maturities[i] + Nmat-1; |
---|
89 | swapval = B[m] + swaprates[i]*S[m] - 1.0; |
---|
90 | condassign(v,-swapval,v-100.0*swapval); |
---|
91 | } |
---|
92 | |
---|
93 | // apply discount // |
---|
94 | |
---|
95 | for (n=0; n<Nmat; n++) |
---|
96 | v = v/(1.0+delta*L[n]); |
---|
97 | delete[](B); |
---|
98 | delete[](S); |
---|
99 | } |
---|
100 | |
---|
101 | /* -------------------------------------------------------- */ |
---|
102 | |
---|
103 | int main(){ |
---|
104 | |
---|
105 | // constants for all threads |
---|
106 | |
---|
107 | // LIBOR interval // |
---|
108 | double delta = 0.25; |
---|
109 | // data for swaption portfolio // |
---|
110 | int Nopt = 15; |
---|
111 | int maturities[] = {4,4,4,8,8,8,20,20,20,28,28,28,40,40,40}; |
---|
112 | double swaprates[] = {.045,.05,.055,.045,.05,.055,.045,.05, |
---|
113 | .055,.045,.05,.055,.045,.05,.055 }; |
---|
114 | |
---|
115 | int i, j, N, Nmat, npath, nthreads, slot; |
---|
116 | double vtot, *v, *lambda, **z, **grad, *gradtot, **xp; |
---|
117 | |
---|
118 | Nmat = 40; |
---|
119 | N = Nmat+40; |
---|
120 | npath = 30; |
---|
121 | nthreads = 5; |
---|
122 | slot = npath/nthreads; |
---|
123 | |
---|
124 | lambda = new double[N]; |
---|
125 | v = new double[npath]; |
---|
126 | gradtot = new double[N]; |
---|
127 | z = new double*[npath]; |
---|
128 | grad = new double*[npath]; |
---|
129 | xp = new double*[npath]; |
---|
130 | for (i=0;i<npath;i++) |
---|
131 | { |
---|
132 | z[i] = new double[N]; |
---|
133 | grad[i] = new double[N+Nmat]; |
---|
134 | xp[i] = new double[N+Nmat]; |
---|
135 | } |
---|
136 | |
---|
137 | for (i=0;i<N;i++) |
---|
138 | { |
---|
139 | gradtot[i] = 0.0; |
---|
140 | lambda[i] = 0.2; |
---|
141 | } |
---|
142 | |
---|
143 | for (j=0; j<npath; j++) |
---|
144 | { |
---|
145 | v[j] = 0; |
---|
146 | for (i=0; i<N; i++) |
---|
147 | xp[j][i]= 0.05; |
---|
148 | for (i=0; i<Nmat; i++) |
---|
149 | { |
---|
150 | z[j][i] = 0.3; |
---|
151 | xp[j][N+i]= 0.3; |
---|
152 | } |
---|
153 | } |
---|
154 | |
---|
155 | omp_set_num_threads(nthreads); |
---|
156 | |
---|
157 | //----------------------------------------------------------// |
---|
158 | // // |
---|
159 | // do a full path + portfolio sensitivity check // |
---|
160 | // // |
---|
161 | // A real application would generate a different random // |
---|
162 | // vector z for each path but here we set one and reuse it // |
---|
163 | // // |
---|
164 | //----------------------------------------------------------// |
---|
165 | |
---|
166 | #pragma omp parallel firstprivate(ADOLC_OpenMP_Handler) |
---|
167 | { |
---|
168 | // different paths for each thread |
---|
169 | int index = omp_get_thread_num(); |
---|
170 | int l,k; |
---|
171 | int rv = 0; |
---|
172 | |
---|
173 | adouble *La, va, *za; |
---|
174 | |
---|
175 | La = new adouble[N]; |
---|
176 | za = new adouble[Nmat]; |
---|
177 | |
---|
178 | int init = index*slot; |
---|
179 | |
---|
180 | trace_on(init); |
---|
181 | for(k=0;k<N;k++) |
---|
182 | La[k] <<= 0.050000; |
---|
183 | for(j=0;j<Nmat;j++) |
---|
184 | za[j] <<= z[init][j]; |
---|
185 | |
---|
186 | path_calc(N,Nmat,delta,La,lambda,za); |
---|
187 | portfolio(N,Nmat,delta,Nopt,maturities,swaprates,La,va); |
---|
188 | |
---|
189 | va >>= v[init]; |
---|
190 | trace_off(1); |
---|
191 | |
---|
192 | for(l=init;l<init+slot;l++) |
---|
193 | { |
---|
194 | rv = gradient(init,N+Nmat,xp[l],grad[l]); |
---|
195 | } |
---|
196 | |
---|
197 | #pragma omp barrier |
---|
198 | if (index==0) |
---|
199 | { |
---|
200 | vtot = 0; |
---|
201 | for (l=0; l<npath; l++) |
---|
202 | { |
---|
203 | vtot += v[l]; |
---|
204 | for(k=0;k<N;k++) |
---|
205 | gradtot[k] += grad[l][k]; |
---|
206 | } |
---|
207 | vtot = vtot/npath; |
---|
208 | for(k=0;k<N;k++) |
---|
209 | gradtot[k] /= npath; |
---|
210 | |
---|
211 | printf("Gradient: \n"); |
---|
212 | for(i=0;i<N;i++) |
---|
213 | printf(" %f \n",gradtot[i]); |
---|
214 | } |
---|
215 | } |
---|
216 | |
---|
217 | |
---|
218 | return 0; |
---|
219 | } |
---|
220 | |
---|