Changeset 143
- Timestamp:
- 04/08/08 21:42:51 (8 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc/model.c (modified) (2 diffs)
- projects/AsynCluster/trunk/svpmc/model.py (modified) (5 diffs)
- projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
projects/AsynCluster/trunk/svpmc/model.c
r141 r143 26 26 // y Observation values, [pn] array 27 27 // a Drift, [p] vector 28 // b Inverted lag-1 cross-correlation [pp] array28 // b Lag-1 cross-correlation [pp] array 29 29 30 30 int i, j, k; … … 83 83 84 84 85 // Model.simulateVolatilities 86 // 87 // Supplied variables 88 // ---------------------------------------------------------------------------- 89 // h Volatility values (initially empty), [pn] array 90 // v Volatility shock values, [pn] array 91 // d Volatility offset, [p] vector 92 // e Lag-1 cross-correlation [pp] array 85 93 94 int i, j, k; 95 double sum; 86 96 97 // The first volatility value for each time series is just the first volatility 98 // shock plus the offset 99 for(i=0; i<Nv[0]; i++) 100 H2(i,0) = V2(i,0) + D1(i); 101 // For each innovation after the first... 102 for(j=1; j<Nv[1]; j++) { 103 // For each time series... 104 for(i=0; i<Nv[0]; i++) { 105 // The offset plus the current shock... 106 sum = D1(i) + V2(i,j); 107 // ...plus the dot product for the VAR(1) term... 108 for(k=0; k<Nv[0]; k++) 109 sum += E2(i,k) * H2(k,j-1); 110 // ...is the modeled output 111 H2(i,j) = sum; 112 } 113 } projects/AsynCluster/trunk/svpmc/model.py
r141 r143 31 31 from twisted_goodies.pybywire.pack import Packer, Unpacker 32 32 33 import weave33 import sample, weave 34 34 35 35 … … 249 249 class VAR(weave.Weaver): 250 250 """ 251 Construct me with a sequence of one or more time-series objects containing252 observation data.251 Construct me with a 2-D array of time series values, one time series per 252 row. 253 253 254 254 Then you can call me with a C{p x 1} drift vector I{a} and a C{p x p} … … 262 262 263 263 """ 264 def __init__(self, tsList): 265 p = len(tsList) 266 tsObj = tsList[0] 267 self.y = s.empty((p, len(tsObj))) 268 for k in xrange(len(tsList)): 269 if k > 0: 270 # No need to intersect first object with itself 271 tsObj = tsObj.intersect(tsList[k]) 272 self.y[k,:] = tsObj() 273 264 def __init__(self, y): 265 self.y = y 266 274 267 def reverse(self, a, b): 275 268 return self.c( … … 281 274 282 275 283 class Model(Parameterized ):276 class Model(Parameterized, weave.Weaver): 284 277 """ 285 278 I am a sophisticated econometric model wherein the residuals of a 286 279 linear-drift, VAR(1) process follow a multivariate stochastic volatility 287 process with correlated volatility and return shocks. 288 280 process with correlated volatility (v) and return (w) shocks. 281 282 @ivar tsList: A sequence of one or more time-series objects containing 283 observation data. 284 289 285 """ 290 286 keyAttrs = {'tsList':None} 291 paramNames = ['a', 'b', 'c', 'd', 'e', 'f'] 292 293 def _get_residualizer(self): 294 if not hasattr(self, '_residualizer'): 295 self._residualizer = Residualizer(self.tsList) 296 return self._residualizer 297 residualizer = property(_get_residualizer) 298 299 def dataToResiduals(self, a, b): 300 """ 301 """ 302 r = self.y 303 weave.inline( 304 """ 305 306 """, 307 ['r', 'a', 'b'], 308 extra_compile_args=['-w'], support_code=self.support) 309 310 def simulateVolatilities(self, svParams): 287 paramNames = [ 288 # Approximately p*(3+2.5) parameters total, e.g., 125 for p=5 289 #------------------------------------------------------------ 290 'a', # Drift (p x 1) 291 292 'b', # Lag-1 cross-correlation for VAR of returns (p x p) 293 294 'c', # Volatility/innovation shock correlations (p x 1) 295 296 'd', # Volatility offset (p x 1) 297 298 'e', # Lag-1 cross-correlations for VAR of volatilities 299 # (p x p) 300 301 'f', # Volatility shock concurrent cross-correlations, 302 # upper triangular of (p x p) 303 ] 304 305 306 #--- Properties ----------------------------------------------------------- 307 308 def _get_p(self): 309 return len(self.tsList) 310 p = property(_get_p) 311 312 def _get_y(self): 313 if 'y' not in self.cache: 314 tsObj = self.tsList[0] 315 y = s.empty((self.p, len(tsObj))) 316 for k in xrange(self.p): 317 if k > 0: 318 # No need to intersect first object with itself 319 tsObj = tsObj.intersect(self.tsList[k]) 320 y[k,:] = tsObj() 321 self.cache['y'] = y 322 return self.cache['y'] 323 y = property(_get_y) 324 325 def _get_n(self): 326 return self.y.shape[1] 327 n = property(_get_n) 328 329 def _get_var(self): 330 if 'var' not in self.cache: 331 self.cache['var'] = VAR(self.tsList) 332 return self.cache['var'] 333 var = property(_get_var) 334 335 def _get_nw(self): 336 if 'nw' not in self.cache: 337 self.cache['nw'] = sample.NormalWalk(self.p, self.n) 338 return self.cache['nw'] 339 nw = property(_get_nw) 340 341 342 #--- Methods -------------------------------------------------------------- 343 344 def simulateVolatilities(self, v): 311 345 """ 312 346 """ 313 347 pass 314 348 315 def rDistribution(self, N): 316 """ 317 Given my current normalization and distribution parameters, returns a 318 1-d array of I{N} random variates from my probability distribution. 319 320 The distribution object caches its result, and repeated calls with the 321 exact same parameters and requested number of samples N will yield the 322 same set of random variates. The reason is to correlate results between 323 other models using the same underlying random variable(s) and 324 distribution(s). 325 """ 326 values = self.distribution.rvs(N) 327 return self.denormalize(values) 328 329 def normalize(self, values=None): 330 """ 331 Override this to do more sophisticated normalization than mere drift 332 removal. 333 334 If no values are supplied, my current data set is used. 335 """ 336 if values is None: 337 values = self.data() 338 if len(self.normParams): 339 return values - sum(self.normParams) 340 return values 341 342 def pDistribution(self, X, N, M=1): 349 350 def pDistribution(self, paramSequence, N, M=1): 343 351 """ 344 352 After application of the normalization and distribution parameters from 345 the supplied parameter vector I{X}, returns a 1d array of I{N} points346 over the range of my normalized data and a 1d array of probability347 densitiesfor each.353 the supplied parameter sequence, returns a 1d array of I{N} points over 354 the range of my normalized data and a 1d array of probability densities 355 for each. 348 356 349 357 You can compute the probability densities over a multiple of the data 350 358 range by setting the keyword I{M} to something greater than 1. 351 359 """ 352 self.setParamVector(X) 353 values = self.data() 354 Y = s.linspace(M*values.min(), M*values.max(), N) 355 Yn = self.normalize(Y) 356 return Y, self.distribution.pdf(Yn) 360 # TODO 357 361 358 362 def likelihood(self, paramSequence): 359 363 """ 360 364 Returns the log-likelihood of my normalized data given my distribution, 361 after application of the normalization and distribution parameters from362 the supplied parameter vector I{X}.365 after application of the parameters from the supplied parameter 366 sequence. 363 367 364 368 This method returns the likelihood of the data given the parameters. It … … 369 373 to this method will be unnecessary. 370 374 """ 371 self.setParamVector(X) 372 Y = self.normalize() 373 return s.log(self.distribution.pdf(Y)).sum() 375 # TODO projects/AsynCluster/trunk/svpmc/test/test_model.py
r141 r143 161 161 class Test_VAR(util.TestCase): 162 162 def setUp(self): 163 self.tsList = [163 tsList = [ 164 164 TimeSeries('prices', PRICES), TimeSeries('earnings', EARNINGS)] 165 self.var = model.VAR(self.tsList) 165 self.y = s.row_stack([ 166 tsList[k].intersect(tsList[1-k])() for k in (0,1)]) 167 self.var = model.VAR(self.y) 166 168 167 169 def _plot(self, *vars):
