Changeset 159

Show
Ignore:
Timestamp:
04/25/08 23:28:27 (7 months ago)
Author:
edsuom
Message:

Innovation & volatility shocks can have non-unity variances; working on integrating PMC into the new scheme of paramContainers

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r158 r159  
    1717 
    1818""" 
    19 The PMC runner 
     19The stochastic volatility model and its manager 
    2020""" 
    2121 
    …… 
    4949 
    5050 
     51class SpecParser(object): 
     52    """ 
     53    """ 
     54    def __init__(self, projectSpec): 
     55        self.spec = projectSpec 
     56 
     57    def priorContainer(self): 
     58        """ 
     59        """ 
     60 
     61    def timeSeries(self): 
     62        """ 
     63        """ 
     64 
     65 
    5166class ModelManager(object): 
    5267    """ 
    …… 
    7590    ########################################################################### 
    7691    """ 
    77     def __init__(self, modelObj): 
    78         self.modelObj = modelObj 
     92    def __init__(self, projectSpec): 
     93        sp = SpecParser(projectSpec) 
     94        self.priors = sp.priorContainer() 
     95        self.modelObj = Model(tsList=sp.timeSeries()) 
    7996        self.queue = NullQueue(1) 
    8097        reactor.addSystemEventTrigger( 
    …… 
    109126        d.addErrback(self._oops) 
    110127        return d 
    111  
     128     
    112129    def likelihood(self, paramContainer, wiggle, remote=False, local=False): 
    113130        """ 
    …… 
    119136        undergoing some nested MCMC. 
    120137 
    121         Returns a deferred that fires (with C{None}) when the likelihood 
    122         computation is done and the paramContainer has been updated. 
     138        Returns a deferred that fires with the log-likelihood when the 
     139        likelihood computation is done and the paramContainer has been updated. 
    123140        """ 
    124141        def gotPackedLikelihood(result): 
    …… 
    128145                paramContainer.pLogLikelihood = -s.inf 
    129146            up = pack.Unpacker(result) 
    130             paramContainer.L = up() 
     147            L = paramContainer.L = up() 
    131148            paramContainer.z = up() 
     149            return L 
    132150 
    133151        def gotLikelihood(result): 
    134             paramContainer.L = result[0] 
     152            L = paramContainer.L = result[0] 
    135153            paramContainer.z = result[1] 
     154            return L 
    136155         
    137156        if (remote or self.remoteMode) and not local: 
    …… 
    226245        if 'rv' not in self.cache: 
    227246            # Concurrent cross-correlations between volatility shocks 
    228             self.cache['rv'] = s.matrix( 
    229                 linalg.cholesky(self.covarMatrix(self.f), lower=True)) 
     247            self.cache['rv'] = s.matrix(linalg.cholesky( 
     248                self.covarMatrix(self.fs, self.fr), lower=True)) 
    230249        return self.cache['rv'] 
    231250    rv = property(_get_rv) 
    …… 
    234253    #--- Support -------------------------------------------------------------- 
    235254 
    236     def covarMatrix(self, correlations): 
    237         """ 
    238         Returns a covariance matrix from the vector or sequence of 
    239         I{correlations} supplied in the form used by L{correlate}. 
    240  
    241         Assumes that the independent components have unit variance. 
     255    def covarMatrix(self, cs, cr): 
     256        """ 
     257        Returns a covariance matrix from the vector or sequence of I{cs} 
     258        standard deviations and I{cr} cross-correlations. 
     259 
     260        The I{cs} argument is in the form s0, s1,..., sp. The I{cr} argument is 
     261        in the form r00, r01, r02,..., r0p, r11, r12,..., r1p,..., rpp 
    242262        """ 
    243263        i = 0 
    244264        p = self.p 
    245         cv = s.ones((p, p)) 
    246         for j in xrange(p - 1): 
     265        cv = s.zeros((p, p)) 
     266        for j in xrange(p): 
     267            cv[j,j] = cs[j]**2 
    247268            for k in xrange(j+1, p): 
    248                 cv[j, k] = cv[k, j] = correlations[i] 
     269                cv[j, k] = cv[k, j] = cr[i] * cs[j] * cs[k] 
    249270                i += 1 
    250271        return cv 
    251  
     272     
    252273    def decaySamples(self, tol): 
    253274        """ 
    …… 
    323344            # Inverse of concurrent cross-correlations between innovation shocks 
    324345            ri = linalg.inv( 
    325                 linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
     346                linalg.cholesky( 
     347                self.covarMatrix(self.cs, self.cr), lower=True)) 
    326348            self.cache['lv_work'] = y, ri 
    327349        # Run the C code where the heavy lifting is done 

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

    r158 r159  
    3838    'b',    # Lag-1 cross-correlation for VAR of returns (p x p) 
    3939 
    40     'c',    # Innovation shock concurrent correlations (vector with 
     40    'cs',   # Innovation shock standard deviations (vector with s0, s1,..., sp) 
     41 
     42    'cr',   # Innovation shock cross-correlation coefficients (vector with 
    4143            # r01, r02,..., r0p, r12,..., r1p,...) 
    4244     
    …… 
    4547    'e',    # Lag-1 cross-correlations for VAR of volatilities 
    4648            # (p x p) 
    47      
    48     'f',    # Volatility shock concurrent correlations (vector with 
    49             # r01, r02,..., r0p, r12,..., r1p,...) 
     49 
     50    'fs',   # Volatility shock standard deviations (vector with s0, s1,..., sp) 
     51     
     52    'fr',    # Volatility shock coefficients (vector with 
     53             # r01, r02,..., r0p, r12,..., r1p,...) 
    5054 
    5155    'g',    # Volatility/innovation shock correlations (p x 1) 
    …… 
    151155class PriorContainer(object): 
    152156    """ 
    153     TODO 
     157    I am a container for array-like instances of L{FlexArray} that contain 
     158    L{Prior} objects for the array elements of my model's parameters. 
    154159    """ 
    155160    paramNames = PARAM_NAMES 
    156161 
    157     def __init__(self, specFilePath): 
    158         fh = open(specFilePath) 
    159         self.spec = fh.read() 
    160         fh.close() 
    161      
     162    def __init__(self, paramSpec): 
     163        # TODO 
     164        pass 
     165 
    162166 
    163167class Prior(Parameterized): 
    …… 
    238242class ParameterContainer(Parameterized): 
    239243    """ 
     244    I am a container for arrays of model parameters. 
    240245    """ 
    241246    keyAttrs = {'z':[], 'L':None} 

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

    r158 r159  
    2121 
    2222import scipy as s 
    23 from scipy import stats, random, linalg 
    24 from scipy.signal import signaltools 
     23from scipy import stats, random, linalg, signal 
    2524 
    2625from twisted.internet import reactor, defer 
    …… 
    166165                "'Uncorrelated' hypothesis not disproved, p=%f" % corrInfo[1]) 
    167166     
     167    def _check_psd(self, x, *coeffs): 
     168        n = x.shape[1] 
     169        fig = self.fig 
     170        for k, coeff in enumerate(coeffs): 
     171            spNumber = 100*len(coeffs) + 21 + 2*k 
     172            sp = fig.add_subplot(spNumber) 
     173            pxx, fxx = sp.psd(x[k,n/2:], NFFT=32) 
     174            hxx = s.absolute( 
     175                signal.freqz([1], [1, -coeff], s.pi*fxx)[1])**2 
     176            sp.plot(fxx, 10*s.log10(hxx)) 
     177            err = (s.absolute(pxx - hxx) / hxx).max() 
     178            sp.set_title( 
     179                "AR coeff = %4.2f, rolloff err = %3.1f " % (coeff, err)) 
     180            sp = fig.add_subplot(spNumber+1) 
     181            sp.plot(x[k,:]) 
     182            self.failUnless(err < 2.5, "Unexpected PSD rollof") 
     183 
    168184 
    169185class Test_VAR(Multivariate_BaseCase): 
    …… 
    192208        self._check_uncorr(y[0,:], y[1,:]) 
    193209        # IIR filtering 
    194         fig = self.fig 
    195         for k, bk in enumerate((b1, b2)): 
    196             sp = fig.add_subplot(221 + 2*k) 
    197             pxx = sp.psd(y[k,n/2:], NFFT=32)[0] 
    198             ratio = ( pxx[-1]/pxx[0] ) / ( 2*((1-bk)/(1+bk))**2 ) 
    199             self.failUnless(ratio > 0.4, "Too little rolloff in PSD") 
    200             self.failUnless(ratio < 1.5, "Too much rolloff in PSD") 
    201             sp = fig.add_subplot(221 + 2*k + 1) 
    202             sp.plot(y[k,:]) 
     210        self._check_psd(y, b1, b2) 
    203211 
    204212    def test_forward_xcorr(self): 
    …… 
    247255 
    248256class Test_Model_misc(util.TestCase): 
     257    def test_covarMatrix(self): 
     258        def tryCoeffs(cs, cr, x): 
     259            x = s.array(x) 
     260            modelObj.tsList = [None] * x.shape[0] 
     261            y = modelObj.covarMatrix(cs, cr) 
     262            self.failUnlessEqual(y.shape, x.shape) 
     263            self.failUnless( 
     264                s.absolute(y-x).max() < 1E-8, 
     265                "Generated array \n%s\n != \n%s" % (y, x)) 
     266         
     267        modelObj = model.Model() 
     268        tryCoeffs([1.0], [], [[1.0]]) 
     269        tryCoeffs([1.0, 1.0], [0.2], [[1.0, 0.2], [0.2, 1.0]]) 
     270        tryCoeffs([2, 3], [0], [[4.0, 0.0], [0.0, 9.0]]) 
     271        tryCoeffs([3, 4], [0.5], [[9.0, 6.0], [6.0, 16.0]]) 
     272        tryCoeffs( 
     273            [3, 4, 5], [0.1, 0.2, 0.3], 
     274            [[9.0, 1.2, 3.0], [1.2, 16.0, 6.0], [3.0, 6.0, 25.0]]) 
     275 
    249276    def test_logUniform(self): 
    250277        N = 35000 
    …… 
    320347 
    321348    def _draw_h(self): 
    322         self.model.c = s.array([0.0]) 
    323         self.model.g = s.array([0.0, 0.0]) 
     349        self.model.cs = s.array([1.0, 1.0]) 
     350        self.model.cr = [0.0] 
     351        self.model.g = s.array([0.0]) 
    324352        z = s.randn(self.model.p, self.model.n) 
    325         rv = s.matrix( 
    326             linalg.cholesky(self.model.covarMatrix(self.model.f), lower=True)) 
     353        rv = s.matrix(linalg.cholesky( 
     354            self.model.covarMatrix(self.model.fs, self.model.fr), lower=True)) 
    327355        v = s.array(rv * z) 
    328356        return self.model.logVolatilities( 
    …… 
    333361        self.model.tsList = [self.model.tsList[0]] 
    334362        n = self.model.n 
    335         self.model.c = [] 
     363        self.model.cs = [1.0] 
     364        self.model.cr = [] 
    336365        self.model.d = s.array([0.0]) 
    337366        self.model.e = s.array([[e]]) 
    338         self.model.f = [] 
     367        self.model.fs = [1.0] 
     368        self.model.fr = [] 
    339369        self.model.g = s.array([0.0]) 
    340370        z = s.randn(1, self.model.n) 
    …… 
    342372            z, z.copy(), s.zeros_like(z), s.array([0.0]))[2] 
    343373        # IIR filtering 
    344         fig = self.fig 
    345         sp = fig.add_subplot(211) 
    346         pxx = sp.psd(h[0,n/2:], NFFT=32)[0] 
    347         ratio = ( pxx[-1]/pxx[0] ) / ( 2*((1-e)/(1+e))**2 ) 
    348         self.failUnless(ratio > 0.4, "Too little rolloff in PSD") 
    349         self.failUnless(ratio < 1.5, "Too much rolloff in PSD") 
    350         sp = fig.add_subplot(212) 
    351         sp.plot(h[0,:]) 
     374        self._check_psd(h, e) 
    352375 
    353376    def test_indep_offset(self): 
    354377        n = self.model.n 
    355378        self.model.d = s.array([1.0, 1.0]) 
    356         self.model.f = s.array([0.0]) 
     379        self.model.fs = [1.0, 1.0] 
     380        self.model.fr = [0.0] 
    357381        for j in xrange(10): 
    358382            e1, e2 = 0.8*s.rand(2) + 0.1 
    …… 
    369393        self.model.d = s.array([0.0, 0.0]) 
    370394        self.model.e = s.array([[e1, 0.0], [0.0, e2]]) 
    371         self.model.f = [0.0] 
     395        self.model.fs = [1.0, 1.0] 
     396        self.model.fr = [0.0] 
    372397        h = self._draw_h() 
    373398        # Uncorrelated 
    374399        self._check_uncorr(h[0,:], h[1,:]) 
    375400        # IIR filtering 
    376         fig = self.fig 
    377         for k, ek in enumerate((e1, e2)): 
    378             sp = fig.add_subplot(221 + 2*k) 
    379             pxx = sp.psd(h[k,n/2:], NFFT=32)[0] 
    380             ratio = ( pxx[-1]/pxx[0] ) / ( 2*((1-ek)/(1+ek))**2 ) 
    381             self.failUnless(ratio > 0.4, "Too little rolloff in PSD") 
    382             self.failUnless(ratio < 1.5, "Too much rolloff in PSD") 
    383             sp = fig.add_subplot(221 + 2*k + 1) 
    384             sp.plot(h[k,:]) 
     401        self._check_psd(h, e1, e2) 
    385402 
    386403    def _check_xcorr(self, dPeak): 
    …… 
    403420        self.model.d = s.array([0.0, 0.0]) 
    404421        self.model.e = s.array([[0.0, 0.9], [0.0, 0.0]]) 
    405         self.model.f = [0.0] 
     422        self.model.fs = [1.0, 1.0] 
     423        self.model.fr = [0.0] 
    406424        self._check_xcorr(1) 
    407425 
    …… 
    409427        self.model.d = s.array([0.0, 0.0]) 
    410428        self.model.e = s.array([[0.0, 0.0], [0.0, 0.0]]) 
    411         self.model.f = [0.5] 
     429        self.model.fs = [1.0, 1.0] 
     430        self.model.fr = [0.5] 
    412431        self._check_xcorr(0) 
    413432 
    …… 
    475494                # Generate some multivariate normal test data having the 
    476495                # generated correlations 
    477                 r = self.model.covarMatrix(c) 
     496                r = self.model.covarMatrix([1.0], c) 
    478497                if VERBOSE: 
    479498                    print "\n\n", p 
    …… 
    502521        'a': [0], 
    503522        'b': [[0]], 
    504         'c': [], 
     523        'cs': [1.0], 
     524        'cr': [], 
    505525        'd': [0], 
    506526        'e': [[0]], 
    507         'f': [], 
     527        'fs': [1.0], 
     528        'fr': [], 
    508529        'g': [0], 
    509530        } 
    …… 
    533554 
    534555class Test_Model_likelihood(Model_BC_Mixin, Multivariate_BaseCase): 
    535     N = 3000 
     556    N = 1000 
    536557    corr = 0.0 
    537558    params = { 
    538559        'a': [0], 
    539560        'b': [[0]], 
    540         'c': [], 
    541         'd': [0.0], 
    542         'e': [[0.95]], 
    543         'f': [], 
     561        'cs': [1.0], 
     562        'cr': [], 
     563        'd': [0], 
     564        'e': [[0]], 
     565        'fs': [1.0], 
     566        'fr': [], 
    544567        'g': [corr], 
    545568        } 
    546569 
    547     class Mock_ParameterContainer(object): 
    548         latent = [] 
    549         def paramSequence(self): 
    550             return [self.params[x] for x in "abcdefg"] 
    551  
    552     @classmethod 
    553     def pcFactory(cls): 
    554         pc = cls.Mock_ParameterContainer() 
    555         pc.params = cls.params 
    556         return pc 
    557      
    558570    def _setupModel(self, **kw): 
    559571        # Any special model parameters 
    …… 
    565577            [[1.0, self.model.g[0]], [self.model.g[0], 1.0]], 
    566578            self.N).transpose() 
    567         self.h = signaltools.lfilter( 
     579        self.h = signal.lfilter( 
    568580            [1], [1.0, -self.model.e[0][0]], self.iv[1,:]) 
    569581        x = s.exp(0.5 * self.h) * self.iv[0,:] 
    …… 
    596608        self._setupModel(e=[[0.95]]) 
    597609        # Sigma=0.5 seems to give us the supposedly ideal ~44% acceptance rate 
    598         self._runHG(50, 0.5) 
     610        self._runHG(20, 0.5) 
    599611         
    600612    def test_hybridGibbs_lowCorrelation(self): 
    601613        self._setupModel(e=[[0.4]]) 
    602         self._runHG(50, 0.5) 
     614        self._runHG(20, 0.5) 
    603615 
    604616    def test_hybridGibbs_noCorrelation(self): 
    605617        self._setupModel(e=[[0.0]]) 
    606         self._runHG(50, 0.5) 
     618        self._runHG(20, 0.5) 

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

    r132 r159  
    4141 
    4242 
    43 class Mock_ModelCaller(util.Mock): 
    44     def __init__(self, modelMap): 
    45         self.modelMap = modelMap 
     43class Mock_ModelManager(util.Mock): 
     44    def __init__(self, modelObj): 
     45        self.modelObj = modelObj 
    4646        self.threadQueue = Mock_ThreadQueue() 
    4747        self.resetCalls() 
    …… 
    5252        self.stepCounter = 0 
    5353         
    54     def call(self, ID, methodName, *args, **kw): 
    55         method = getattr(self.modelMap[ID], methodName) 
    56         return self.threadQueue.call(method, *args, **kw) 
    57  
    58     def likelihood(self, IDs, paramVectors, **kw): 
     54    def likelihood(self, paramContainer, wiggle, remote=False, local=False): 
    5955        def done(results): 
    6056            if self.callsPending: 
    6157                self.stepCounter += 1 
    6258            self.callsPending = [] 
    63             return sum(results) 
    64  
     59            paramContainer.L = results[0] 
     60            paramContainer.z = results[1] 
     61            return results[0] 
     62         
    6563        self.callCounter += 1 
    6664        self.callsPending.append(self.callCounter) 
    67         dList = [ 
    68             self.threadQueue.call(self.modelMap[ID].likelihood, paramVectors[k]) 
    69             for k, ID in enumerate(IDs)] 
    70         return defer.gatherResults(dList).addCallback(done) 
    71  
    72  
    73 class Mock_ProjectManager(util.Mock): 
    74     expressions = {} 
    75     paramNames, priors, models = util.gaussianModelFactory() 
     65        d = self.threadQueue.call(self.modelObj.likelihood, paramContainer) 
     66        d.addCallback(done) 
     67        return d 
    7668 
    7769 
    …… 
    109101     
    110102    def setUp(self): 
    111         projectManager = Mock_ProjectManager() 
    112103        self.mgr = model.ModelManager(projectManager) 
    113104        self.mgr.caller = Mock_ModelCaller(self.mgr.modelMap) 
    …… 
    130121        ri.rc('plot', mo) 
    131122        return util.deferToLater(delay=delay) 
    132  
    133  
    134 class Test_MCMC(BaseTC): 
    135     def test_getXL(self): 
    136         def gotPopulation(XL): 
    137             self.failUnlessEqual(len(XL), 1000) 
    138             X, Y, Z = XL[:,0], XL[:,1], XL[:,2] 
    139             I = s.isfinite(Z).nonzero() 
    140             ax = p3.Axes3D(self.fig) 
    141             ax.scatter3d(X[I], Y[I], Z[I]) 
    142  
    143         return self.mc.getXL(1000).addCallback(gotPopulation) 
    144  
    145     @defer.deferredGenerator 
    146     def test_metropolisHastings(self): 
    147         N_iter = 500 
    148         self.mc.XL = s.array([[0.1, 0.9, 0]]) 
    149         acceptances = 0 
    150         for i in xrange(N_iter): 
    151             wfd = defer.waitForDeferred(self.mc.metropolisHastings(0)) 
    152             yield wfd 
    153             wasAccepted = wfd.getResult() 
    154             acceptances += wasAccepted 
    155         print "\nAccepted %d of %d proposals" % (acceptances, N_iter) 
    156         for func, text in ((lambda x: x > 0.3*N_iter, "low"), 
    157                            (lambda x: x < 0.6*N_iter, "high")): 
    158             self.failUnless( 
    159                 func(acceptances), 
    160                 "Acceptance rate of %d%% is too %s" \ 
    161                 % (100.0*acceptances/N_iter, text)) 
    162  
    163     def test_run_convergence(self): 
    164         self.mc.N_iter = 1000 
    165         label = "Single chain, %d iterations" % self.mc.N_iter 
    166         self.mc.subscribe(self.consumer) 
    167         d = self.mc.run(1) 
    168         d.addCallback(self._analyzeSamples) 
    169         return d 
    170  
    171  
    172 class Test_DE_MCMC(BaseTC): 
    173     def test_getXds(self): 
    174         def gotPopulation(XL): 
    175             self.mc.XL = XL 
    176             Xds = self.mc.getXds() 
    177             print Xds 
    178  
    179         self.mc.N_chains = 10 
    180         return self.mc.getXL(self.mc.N_chains).addCallback(gotPopulation) 
    181  
    182     def test_run_concurrence(self): 
    183         self.mc.N_iter, N_chains = 10, 100 
    184         self.mc.subscribe(self.consumer) 
    185         return self.mc.run(N_chains) 
    186  
    187     def test_run_convergence(self): 
    188         self.mc.N_iter, N_chains = 100, 20 
    189         label = "Population of %d members, %d generations" \ 
    190                 % (N_chains, self.mc.N_iter) 
    191         self.mc.subscribe(self.consumer) 
    192         d = self.mc.run(N_chains) 
    193         d.addCallback(self._analyzeSamples) 
    194         return d 
    195123 
    196124