Changeset 147

Show
Ignore:
Timestamp:
04/16/08 16:27:16 (8 months ago)
Author:
edsuom
Message:

Got decorrelation working; weave.Weaver now uses inline as its method name, like weave.inline

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/model.c

    r145 r147  
    112112  } 
    113113} 
     114 
     115 
     116 
     117// Model.decorrelate 
     118// 
     119// Supplied variables 
     120// ---------------------------------------------------------------------------- 
     121// y    2-D array of correlated random variates (to be overwritten) 
     122// x    2-D array of independent random variates 
     123// r    2-D array of inverted correlations between values in each column of x 
     124 
     125int i, j, k; 
     126double sum; 
     127 
     128// For each column of independent variates to be decorrelated... 
     129for(i=0; i<Nx[1]; i++) { 
     130  // For each of the values in the column... 
     131  for(j=0; j<Nx[0]; j++) { 
     132    // The product of the inverted cross-correlation matrix and the correlated 
     133    // variates... 
     134    sum = 0; 
     135    for(k=0; k<Nx[0]; k++) 
     136      sum += R2(j,k) * X2(k,i); 
     137    // ...is the correlated output 
     138    Y2(j, i) = sum; 
     139  } 
     140} 

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

    r146 r147  
    266266     
    267267    def reverse(self, a, b): 
    268         return self.c('v', 'y', 'a', 'b', v=s.empty_like(self.y)) 
     268        return self.inline('v', 'y', 'a', 'b', v=s.empty_like(self.y)) 
    269269 
    270270    def forward(self, v, a, b): 
    271         return self.c('y', 'v', 'a', 'b', y=s.empty_like(v)) 
     271        return self.inline('y', 'v', 'a', 'b', y=s.empty_like(v)) 
    272272 
    273273 
    …… 
    290290        'b',    # Lag-1 cross-correlation for VAR of returns (p x p) 
    291291 
    292         # The attribute name 'c' is used by the Weaver inline method 
     292        'c',    # Innovation shock concurrent correlations (vector with 
     293                # r01, r02,..., r0p, r12,..., r1p,...) 
    293294         
    294295        'd',    # Volatility offset (p x 1) 
    …… 
    300301                # r01, r02,..., r0p, r12,..., r1p,...) 
    301302 
    302         #--- TODO ---------------- 
    303303        'g',    # Volatility/innovation shock correlations (p x 1) 
    304         #------------------------- 
    305304 
    306305        ] 
    …… 
    354353            self.h = s.empty_like(self.nw.x) 
    355354        self.nw.step(wiggle) 
    356         return self.c('h', 'v', 'd', 'e', v=self.nw.correlate(self.f)) 
     355        return self.inline('h', 'v', 'd', 'e', v=self.nw.correlate(self.f)) 
     356 
     357    def decorrelate(self, x, h=None): 
     358        """ 
     359        Call this method with a C{[pn]} array I{x} containing one row of C{n} 
     360        innovations for each of C{p} Wiener processes, with unit variance 
     361        assumed and correlated by a covariance matrix constructed from my 
     362        correlation parameters I{c}. The result will be the decorrelated 
     363        innovations of the Wiener processes. 
     364 
     365        You can supply a C{[pn]} array I{h} of log-volatilities for each of the 
     366        innovations. Then the innovations I{x} are assumed to be scaled by 
     367        their respective volatilities, and will be inversely scaled to meet the 
     368        assumption of unit variance. 
     369        """ 
     370        xs = x.shape 
     371        if xs[0] != self.p: 
     372            raise ValueError("Data array must have %d rows" % self.p) 
     373        hs = s.shape(h) 
     374        if hs: 
     375            if hs != xs: 
     376                raise ValueError("Need one log-volatility per innovation") 
     377            x *= s.exp(-0.5*h) 
     378        if 'decorrelate' not in self.cache: 
     379            self.cache['decorrelate'] = linalg.inv( 
     380                linalg.cholesky(self.nw.covarMatrix(self.c), lower=True)) 
     381        r = self.cache['decorrelate'] 
     382        return self.inline('y', 'x', 'r', y=s.empty_like(x)) 
     383 
     384    def pInnovations(self, y): 
     385        """ 
     386        Call this method with a C{[pn]} array I{y} of uncorrelated 
     387        innovations. 
     388 
     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        """ 
    357393 
    358394 
    359395    #--- The API -------------------------------------------------------------- 
    360396 
    361     def likelihood(self, wiggle, paramSequence): 
    362         """ 
    363         Returns a log-likelihood of my observations for one set of simulated 
    364         volatilies, given the supplied fractional I{wiggle} and the parameter 
    365         arrays in the supplied sequence I{paramSequence}. 
    366  
    367         This method returns the likelihood of the data given the parameters. It 
    368         does not consider the prior probability of the parameters. You have to 
    369         do that yourself, preferably before calling this method. If the prior 
    370         probability of any parameter is zero, the posterior density for that 
    371         parameter will also be zero, and a call to this method will be 
    372         unnecessary. 
     397    def likelihood(self, paramSequence, v, wiggle): 
     398        """ 
     399        Call this method with a sequence of arrays I{paramSequence} defining my 
     400        model parameters and a C{[pn]} array of volatility shocks I{v} to be 
     401        used as a starting point for another simulation draw of volatilities. 
     402 
     403        A new 2-D array of log-volatilities will be drawn from a relatively 
     404        short Metropolis-Hastings random walk simulation of the volatility 
     405        shocks coupled with a VAR(1) from the shocks to the 
     406        log-volatilities. The Metropolis-Hastings proposals are scaled by the 
     407        supplied fractional I{wiggle}. The simulation is governed by the 
     408        log-likelihood of my observation data given the supplied parameters and 
     409        the simulated volatilities. 
     410 
     411        The value of that log-likelihood at the end of the log-volatility 
     412        simulation is returned along with the with the volatility shocks 
     413        underlying the log-volatilities drawn at that point. 
     414 
     415        The returned likelihood does not consider the prior probability of the 
     416        parameters. You have to account for that yourself, preferably before 
     417        calling this method. If the prior probability of any parameter is zero, 
     418        the posterior density for that parameter will also be zero no matter 
     419        what, making any call to this method that uses it a waste of computing 
     420        time. 
    373421        """ 
    374422        self.setParamSequence(paramSequence) 
    375423        # Obtain innovations as residuals of the reversed VAR(1) 
    376424        x = self.var.reverse(self.a, self.b) 
    377         # Draw a new set of simulated volatilities 
    378         h = self.draw_h(wiggle) 
     425        # Simulate log-volatilities for a while 
     426        h = self.draw_h(v, wiggle) 
    379427        # Compute the conditional densities 
    380428         

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

    r146 r147  
    6363        x[:,0] = K*W[I0] 
    6464        a = s.less(x[:,0], 1.0).nonzero()[0] 
    65         x = self.c('x', 'a', 'b', b=s.setdiff1d(s.arange(0, K), a)) 
     65        x = self.inline('x', 'a', 'b', b=s.setdiff1d(s.arange(0, K), a)) 
    6666        V = s.rand(N) 
    6767        RI = s.random.randint(0, K, N) 
    …… 
    8989        A reference to the updated array is returned. 
    9090        """ 
    91         return self.c('x', 'p', 'm', 'wiggle') 
     91        return self.inline('x', 'p', 'm', 'wiggle') 
    9292 
    93     def _covarMatrix(self, correlations): 
     93    def covarMatrix(self, correlations): 
    9494        """ 
    95         Returns a covariance matrix from the supplied vector or sequence of 
    96         I{correlations}, in the form used by L{correlate}. 
     95        Returns a covariance matrix from the vector or sequence of 
     96        I{correlations} supplied in the form used by L{correlate}. 
     97 
     98        Assumes that the independent components have unit variance, as is the 
     99        case with my random walkers C{x}. 
    97100        """ 
    98101        i = 0 
    …… 
    124127            # next time 
    125128            self.r = linalg.cholesky( 
    126                 self._covarMatrix(correlations), lower=True) 
     129                self.covarMatrix(correlations), lower=True) 
    127130            self.y = s.empty_like(self.x) 
    128131        n = self.x.shape[0] 
    …… 
    131134                "Must have a [%d x %d] cross-correlation matrix defined" \ 
    132135                % (n, n)) 
    133         return self.c('y', 'x', 'r', y=s.empty_like(self.x)) 
     136        return self.inline('y', 'x', 'r', y=s.empty_like(self.x)) 
    134137         
    135138         

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

    r146 r147  
    280280        self.model.f = [0.5] 
    281281        self._check_xcorr(0) 
     282 
     283 
     284class Test_Model_decorrelate(util.TestCase): 
     285    def setUp(self): 
     286        self.model = model.Model() 
     287 
     288    def test_basic(self): 
     289        n = 1000 
     290        for p in xrange(2, 6): 
     291            for correlation in s.linspace(0.1, 0.9, 10): 
     292                # Generate some multivariate normal test data 
     293                r = s.maximum(s.identity(p), correlation*s.ones((p,p))) 
     294                x = random.multivariate_normal(s.zeros(p), r, n).transpose() 
     295                # Tweak the model object to make it testable for this p value 
     296                self.model.tsList = [util.TS_LIST[0]] * p 
     297                self.model.c = correlation * s.ones( 
     298                    sum([p-k-1 for k in xrange(p)])) 
     299                # Get the decorrelated values 
     300                y = self.model.decorrelate(x) 
     301                # Check for non-correlation 
     302                for j in xrange(p): 
     303                    for k in xrange(p): 
     304                        if j == k: 
     305                            continue 
     306                        corrInfo = stats.spearmanr(y[j,:], y[k,:]) 
     307                        self.failUnless( 
     308                            abs(corrInfo[0]) < 0.1, 
     309                            "Doesn't appear uncorrelated") 

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

    r145 r147  
    116116    def test_covarMatrix_2x2(self): 
    117117        walker = sample.NormalWalk(2, 2) 
    118         x = walker._covarMatrix([0.5]) 
     118        x = walker.covarMatrix([0.5]) 
    119119        y = s.array([[1, 0.5], [0.5, 1.0]]) 
    120120        self.failUnless(s.equal(x, y).all()) 
    …… 
    122122    def test_covarMatrix_3x3(self): 
    123123        walker = sample.NormalWalk(3, 67) 
    124         x = walker._covarMatrix([0.12, 0.13, 0.23, 0.24]) 
     124        x = walker.covarMatrix([0.12, 0.13, 0.23, 0.24]) 
    125125        y = s.array([ 
    126126            [1.00, 0.12, 0.13], 
    …… 
    131131    def test_covarMatrix_4x4(self): 
    132132        walker = sample.NormalWalk(4, 100) 
    133         x = walker._covarMatrix([0.12, 0.13, 0.14, 0.23, 0.24, 0.34]) 
     133        x = walker.covarMatrix([0.12, 0.13, 0.14, 0.23, 0.24, 0.34]) 
    134134        y = s.array([ 
    135135            [1.00, 0.12, 0.13, 0.14], 

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

    r145 r147  
    5656 
    5757    def __call__(self, x): 
    58         return self.c('x')[0] 
     58        return self.inline('x')[0] 
    5959     
    6060    def foo(self, x): 
    61         return self.c('x', 'y', x=2*s.array([x]))[0] 
     61        return self.inline('x', 'y', x=2*s.array([x]))[0] 
    6262 
    6363 
    …… 
    6666        self.weaver = TestableWeaver() 
    6767 
    68     def test_c_alone(self): 
     68    def test_inline_alone(self): 
    6969        self.weaver._code = {'support':None, 'foo':CODE} 
    7070        self.failUnlessEqual(self.weaver.foo(s.array([1])), 1164) 
    …… 
    7474        self.failUnlessElementsEqual(code.keys(), ['support', 'foo', '__call__']) 
    7575 
    76     def test_c_normal(self): 
     76    def test_inline_normal(self): 
    7777        self.weaver._code = self.weaver._parseCode(CODE) 
    7878        self.failUnlessEqual(self.weaver.foo(10), 625) 
    7979 
    80     def test_c_privateMethod(self): 
     80    def test_inline_privateMethod(self): 
    8181        self.weaver._code = self.weaver._parseCode(CODE) 
    8282        self.failUnlessEqual(self.weaver(s.array([7])), 107) 

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

    r145 r147  
    7676    code = property(_get_code) 
    7777 
    78     def c(self, *args, **kw): 
     78    def inline(self, *args, **kw): 
    7979        """ 
    8080        Call this method to run the inline code for the calling method of the