| 389 | | The result will be the probability density of each innovation, given my |
|---|
| 390 | | volatility/innovation shock correlations C{g} and the volatility shocks |
|---|
| 391 | | currently in my random walker C{nw}. |
|---|
| 392 | | """ |
|---|
| 393 | | |
|---|
| | 394 | The result will be the total log-likelihood of the innovations, given |
|---|
| | 395 | my volatility/innovation shock correlations C{g} and the volatility |
|---|
| | 396 | shocks currently in my random walker C{nw}. A constant term in the |
|---|
| | 397 | log-likelihood is disregarded, C{sqrt(2*pi)}. |
|---|
| | 398 | """ |
|---|
| | 399 | mu = self.g * self.nw.x |
|---|
| | 400 | sigma2 = 1 - self.g**2 |
|---|
| | 401 | logProbs = -(y - mu)**2 / (2*sigma2) - 0.5 * s.log(sigma2) |
|---|
| | 402 | return logProbs.sum() |
|---|
| | 403 | |
|---|
| | 404 | def logUniform(self, N): |
|---|
| | 405 | """ |
|---|
| | 406 | Returns a log-uniform variate taken from a large array that is |
|---|
| | 407 | generated and refreshed periodically. Generating the uniform variates |
|---|
| | 408 | and taking their log in a single large array operation is more |
|---|
| | 409 | efficient than doing so one value at a time. |
|---|
| | 410 | """ |
|---|
| | 411 | logU = getattr(self, '_logU_array', []) |
|---|
| | 412 | if self._logU_index >= len(logU): |
|---|
| | 413 | self._logU_index = -1 |
|---|
| | 414 | # The efficient array operation |
|---|
| | 415 | logU = self._logU_array = s.log(s.rand(10000)) |
|---|
| | 416 | self._logU_index += 1 |
|---|
| | 417 | return logU[self._logU_index] |
|---|
| 426 | | # Simulate log-volatilities for a while |
|---|
| 427 | | h = self.draw_h(v, wiggle) |
|---|
| 428 | | # Compute the conditional densities |
|---|
| 429 | | p = self.pInnovations(self.decorrelate(x, h)) |
|---|
| 430 | | # Return the log-likelihood and the underlying volatility shocks |
|---|
| 431 | | return s.exp(p).sum(), self.nw.x |
|---|
| | 450 | # MCMC simulation of log-volatilities |
|---|
| | 451 | L = -s.inf |
|---|
| | 452 | for j in xrange(self.M): |
|---|
| | 453 | # Draw a proposal array of log-volatilities |
|---|
| | 454 | h = self.draw_h(v, wiggle) |
|---|
| | 455 | # Compute the total log-likelihood conditional on the proposed |
|---|
| | 456 | # log-volatilities |
|---|
| | 457 | Lp = self.likelihoodOfInnovations(self.decorrelate(x, h)) |
|---|
| | 458 | # Metropolis-Hastings accept/reject of the proposal |
|---|
| | 459 | if Lp > L or (Lp-L) > self.logUniform(): |
|---|
| | 460 | L = Lp |
|---|
| | 461 | vShocks = self.nw.x.copy() |
|---|
| | 462 | return L, vShocks |
|---|
