-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathTimeSeries.hpp
427 lines (380 loc) · 6.63 KB
/
TimeSeries.hpp
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
#pragma once
#include <armadillo>
using namespace arma;
namespace TimeSeries {
constexpr double pmSd = 3.0;
constexpr double minPr = 1.0E-4;
struct DiscreteRV {
vec support;
vec probabilities;
};
struct MvNormal {
vec mean;
mat varCovar;
};
/*
* Normal helpers
*/
DiscreteRV
discretizeNormal(const double mean,
const double sigma,
const uword n_elem = 21,
const double pmSigma = pmSd,
bool cutTails = true);
vec discretizeNormal(vec grid,
const double mean,
const double sigma,
bool cutTails = true);
/*
*
* AR(1) process
*
*/
class AR {
public:
inline double
intercept() const
{
return m_intercept;
}
inline double
rho() const
{
return m_rho;
}
inline double
sigma() const
{
return m_sigma;
}
inline double
stationaryMean() const
{
return m_intercept / (1.0 - m_rho);
}
inline double
stationarySigma() const
{
return m_sigma / std::sqrt(1.0 - m_rho * m_rho);
}
inline bool
stationaryQ() const
{
return std::abs(m_rho) < 1.0;
}
AR(double intercept, double rho, double sigma)
: m_intercept(intercept)
, m_rho(rho)
, m_sigma(sigma)
{
}
private:
double m_intercept;
double m_rho;
double m_sigma;
};
/*
*
* VAR(1) process
*
*/
class VAR {
public:
inline uword
size() const
{
return m_size;
}
inline const vec&
intercept() const
{
return m_intercept;
}
inline const mat&
rho() const
{
return m_rho;
}
inline const mat&
sigma() const
{
return m_sigma;
}
const MvNormal conditional(vec& current) const;
vec& stationaryMean();
mat& stationarySigma();
bool stationaryQ();
void print();
VAR(vec intercept, mat rho, mat sigma)
: m_intercept(intercept)
, m_rho(rho)
, m_sigma(sigma)
{
m_size = intercept.n_elem;
}
private:
uword m_size;
vec m_intercept;
mat m_rho;
mat m_sigma;
vec m_statMean;
mat m_statSigma;
bool is_stationary;
void findStationary();
};
/*
*
* Orthogonalized VAR(1) and "rotation" matirx
*
*/
enum OrthoMethod { //
SVD,
Cholesky
};
class OrthogonalizedVAR {
public:
inline const VAR&
getVAR() const
{
return m_ortho;
}
inline const VAR&
getOriginalVAR() const
{
return m_original;
}
inline const mat&
getSupportRotationMatrix() const
{
return LL;
}
OrthogonalizedVAR(VAR original, OrthoMethod method = OrthoMethod::Cholesky)
: m_original(original)
, m_ortho(original)
{
const mat mEye = eye(original.size(), original.size());
mat Atilde = original.intercept();
mat Btilde = original.rho();
mat GG = original.sigma();
vec DD = ones(original.size());
mat diagDD = mEye;
LL = mEye;
switch (method) {
case OrthoMethod::SVD: {
mat VV;
svd(LL, DD, VV, GG);
Atilde = LL.t() * original.intercept();
Btilde = LL.t() * original.rho() * LL;
diagDD = diagmat(DD);
break;
}
case OrthoMethod::Cholesky: {
LL = chol(GG, "lower");
mat invLL = inv(LL);
Atilde = invLL * original.intercept();
Btilde = invLL * original.rho() * LL;
break;
}
default: {
cout << "Unknown method." << endl;
exit(1);
}
}
m_ortho = VAR(Atilde, Btilde, diagDD);
}
private:
VAR m_original;
VAR m_ortho;
mat LL;
};
/*
*
* Markov Chain
*
*/
rowvec
stationaryDistribution(const mat& transitionMatrix);
uvec simulateChain(const mat& transitionMatrix,
const uword initState,
const uword simSz);
class MarkovChain {
public:
inline uword
size() const
{
return m_size;
}
inline rowvec&
support()
{
return m_support;
}
inline const rowvec&
support() const
{
return m_support;
}
inline mat&
transition()
{
return m_tran;
}
inline const mat&
transition() const
{
return m_tran;
}
void save(std::string fname) const;
void print() const;
const rowvec& stationary();
MarkovChain() { }
MarkovChain(rowvec support, mat transition)
: m_support(support)
, m_tran(transition)
{
m_size = support.n_elem;
}
MarkovChain(const AR process,
const uword sz = 21,
const double pmMCsd = pmSd,
bool expSupport = false);
private:
uword m_size;
rowvec m_support;
rowvec m_stationary;
mat m_tran;
};
/*
*
* Discrete approximation of VAR(1)
*
*/
class DiscreteVAR {
public:
inline uword
size() const
{
return m_size;
}
inline uword
flatSize() const
{
return m_flatSize;
}
inline uword
midIx() const
{
return m_midIx;
}
inline const vec
grid(uword varIx) const
{
return m_grids.col(varIx);
}
inline const rowvec
gridValues(uword ix) const
{
return m_grids.row(ix);
}
inline const mat
transition() const
{
return m_tran;
}
void save(std::string fname) const;
void print() const;
DiscreteVAR(const VAR var,
const uword supportSize,
bool trimGrids = true,
OrthoMethod mm = OrthoMethod::Cholesky)
: m_var(var)
{
m_supportSizes.set_size(m_var.size());
m_supportSizes.fill(supportSize);
impl(trimGrids, mm);
}
DiscreteVAR(const VAR var,
const uvec gridSizes,
bool trimGrids = true,
OrthoMethod mm = OrthoMethod::Cholesky)
: m_var(var)
, m_supportSizes(gridSizes)
{
impl(trimGrids, mm);
}
private:
VAR m_var;
uword m_size;
uword m_flatSize;
uword m_midIx;
mat m_grids;
uvec m_supportSizes;
mat m_tran;
void impl(bool trimGrids, OrthoMethod mm);
};
/*
*
* Discrete approximation to VAR(1) with common AR(1) stochastic volatility,
* Basal-Yaron-like.
*
*/
class DiscreteStochVolVAR {
public:
inline uword
size() const
{
return m_size;
}
inline uword
flatSize() const
{
return m_flatSize;
}
inline uword
midIx() const
{
return m_midIx;
}
inline const mat
transition() const
{
return m_tran;
}
inline const vec
grid(uword iIx) const
{
return m_grids.col(iIx);
}
inline const rowvec
gridValues(uword ix) const
{
return m_grids.row(ix);
}
void save(std::string fname) const;
void print() const;
DiscreteStochVolVAR(const VAR var,
const AR vol,
const uvec gridSizes,
const uword volGridSize,
bool trimGrids = true)
: m_var(var)
, m_vol(vol)
, m_supportSizes(gridSizes)
, m_volGridSize(volGridSize)
{
impl(trimGrids);
};
private:
VAR m_var;
AR m_vol;
uword m_size;
uword m_flatSize;
uword m_midIx;
mat m_grids;
mat m_tran;
uvec m_supportSizes;
uword m_volGridSize;
void impl(bool trimGrids);
};
void trimMarkovChain(mat& grids, mat& transition);
} // namespace TimeSeries