Changeset 136

Show
Ignore:
Timestamp:
04/05/08 12:15:17 (8 months ago)
Author:
edsuom
Message:

Working on MCMC simulation of the normals underlying the volatilities

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/sample.py

    r135 r136  
    2222import scipy as s 
    2323from scipy import random, weave 
    24 from scipy.stats.distributions as dists 
     24from scipy.stats import distributions as dists 
    2525 
    2626 
    …… 
    8888class NormalWalk(object): 
    8989    """ 
    90     What is a good approximation to the inverse of the cumulative normal 
    91     distribution function? 
    92      
    93     Mathematically we can write this as:  Find X such that Q(X) = p for 
    94     any 0 < p < 1. (Note: this computes an upper tail probability.) 
    95      
    96     Again, it is not possible to write this as a closed form expression, 
    97     so we resort to approximations. Because of the symmetry of the normal 
    98     distribution, we need only consider 0 < p < 0.5. If you have p > 0.5, 
    99     then apply the algorithm below to q = 1-p, and then negate the value 
    100     for X obtained. (This approximation is also from Abramowitz and 
    101     Stegun.) 
    102      
    103     t = sqrt[ ln(1/p^2) ] 
    104      
    105                  c_0 + c_1*t + c_2*t^2 
    106     X = t -  ------------------------------ 
    107              1 + d_1*t + d_2*t^2 + d_3*t^3 
    108      
    109     c_0 = 2.515517 
    110     c_1 = 0.802853 
    111     c_2 = 0.010328 
    112      
    113     d_1 = 1.432788 
    114     d_2 = 0.189269 
    115     d_3 = 0.001308 
    116      
    117     See Abramowitz and Stegun; Press, et al. 
     90    I run an array of values through a random walk within a zero-mean, unit 
     91    variance normal distribution. Instantiate me with the dimensions of the 
     92    array as one or more arguments. 
    11893    """ 
     94    def __init__(self, *dims): 
     95        self._dims = dims 
     96        self.distObj = dists.norm() 
     97        self.X = s.rand(*dims) 
     98 
    11999    def cinv(self, p): 
    120100        """ 
    …… 
    122102 
    123103        Returns the random variate value(s) occurring with probability less 
    124         than or equal to the specified probabilit(ies) from 0.0 - 1.0. 
     104        than or equal to the specified probabilit(ies) from 0.0 to +1.0. 
    125105        """ 
    126         # TODO: Write C-weave code to do this lightning-fast with the above 
    127         # approximation. 
    128         if not hasattr(self, 'distObj'): 
    129             self.distObj = dists.norm() 
     106        # TODO: Write C-weave code to do this lightning-fast with the following 
     107        # approximation. Because of the symmetry of the normal distribution, we 
     108        # need only consider 0 < p < 0.5. If you have p > 0.5, then apply the 
     109        # algorithm below to q = 1-p, and then negate the value for X obtained. 
     110        # 
     111        # t = sqrt[ ln(1/p^2) ] 
     112        # 
     113        #             c_0 + c_1*t + c_2*t^2 
     114        # X = t -  ------------------------------ 
     115        #          1 + d_1*t + d_2*t^2 + d_3*t^3 
     116        # 
     117        # c_0 = 2.515517 
     118        # c_1 = 0.802853 
     119        # c_2 = 0.010328 
     120        #  
     121        # d_1 = 1.432788 
     122        # d_2 = 0.189269 
     123        # d_3 = 0.001308 
     124        # 
     125        # See Abramowitz and Stegun; Press, et al. 
    130126        return self.distObj.ppf(p) 
     127     
     128    def step(self, W): 
     129        """ 
     130        Call this method with an array I{W} of fractional wiggle values in the 
     131        range 0.0 to +1.0, with one array element for each of my random 
     132        walkers. The elements of my value array will be perturbed accordingly. 
    131133 
    132     def kernel(self, x, sdev): 
     134        Two like-dimensioned arrays are returned, one with the target values 
     135        after the jumps and another containing the probability densities of the 
     136        jumps. 
    133137        """ 
    134         For each current random variate value in I{x}, returns the transition 
    135         kernel... 
    136         """ 
     138        if s.shape(W) != self._dims: 
     139            raise ValueError( 
     140                "Wiggle array must have same dimensions as walker array") 
     141        # TODO: Make work with 2-D array 
     142        P = s.empty_like(W) 
     143        for k, dx in enumerate(0.5*W): 
     144            distObj = dists.truncnorm(self.X[k]-dx, self.X[k]+dx) 
     145            self.X[k] = distObj.rvs(1)[0] 
     146            P[k] = distObj.pdf(self.X[k]) 
     147        return self.X, P 
     148     
     149             
     150         

  • projects/AsynCluster/trunk/svpmc/test/test_sample.py

    r135 r136  
    6262        sp = self.fig.add_subplot(111) 
    6363        sp.hist(Y, bins=20) 
     64 
     65 
     66class Test_NormalWalk(util.TestCase): 
     67    def test_univariate(self): 
     68        walker = sample.NormalWalk(1) 
     69        W = s.array([0.2]) 
     70        X = s.array([walker.step(W)[0][0] for k in xrange(50000)]) 
     71        fig = self.fig 
     72        sp = fig.add_subplot(211) 
     73        sp.plot(X) 
     74        sp = fig.add_subplot(212) 
     75        sp.hist(X, bins=20) 
     76 
     77             
     78