Changeset 141

Show
Ignore:
Timestamp:
04/08/08 03:38:21 (9 months ago)
Author:
edsuom
Message:

Got VAR(1)+drift process working forward & back

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r140 r141  
    2323// Supplied variables 
    2424// ---------------------------------------------------------------------------- 
    25 // r    Residual values (initially empty), [pn] array 
     25// v    Innovation values (initially empty), [pn] array 
    2626// y    Observation values, [pn] array 
    27 // a    Drift vector, [p] vector 
     27// a    Drift, [p] vector 
     28// b    Inverted lag-1 cross-correlation [pp] array 
     29 
     30int i, j, k; 
     31double sum; 
     32 
     33// The first innovation for each time series is just the first observation 
     34// minus the drift 
     35for(i=0; i<Ny[0]; i++) 
     36  V2(i,0) = Y2(i,0) - A1(i); 
     37// For each observation after the first... 
     38for(j=1; j<Ny[1]; j++) { 
     39  // For each time series... 
     40  for(i=0; i<Ny[0]; i++) { 
     41    // The current observation minus the drift... 
     42    sum = Y2(i,j) - A1(i); 
     43    // ...minus the dot product for the VAR(1) term... 
     44    for(k=0; k<Ny[0]; k++) 
     45      sum -= B2(i,k) * Y2(k,j-1); 
     46    // ...is the innovation that yielded this observation 
     47    V2(i,j) = sum; 
     48  } 
     49} 
     50 
     51 
     52 
     53// VAR.forward 
     54//  
     55// Supplied variables 
     56// ---------------------------------------------------------------------------- 
     57// y    Modeled values (initially empty), [pn] array 
     58// v    Innovation values, [pn] array 
     59// a    Drift, [p] vector 
    2860// b    Lag-1 cross-correlation [pp] array 
    2961 
    …… 
    3163double sum; 
    3264 
    33 // For each time series... 
    34 for(i=0; i<Ny[0]; i++) { 
    35   // The model is assumed to fit perfectly until the first observation 
    36   R1(i,0) = 0; 
    37   // For each observation after the first... 
    38   for(j=1; j<Ny[1]; j++) { 
    39     // Start with the inverse drift and the current observation index... 
    40     sum = Y1(i,j) - A1(i); 
    41     // ...then subtract the dot product for the VAR(1) term 
    42     for(k=0; k<Ny[0]; k++) 
    43       sum -= B(i,k) * Y1(k,j-1); 
    44     // The result is the residual for this observation 
    45     R1(i,j) = sum; 
     65// The first modeled value for each time series is just the first innovation 
     66// plus the drift 
     67for(i=0; i<Nv[0]; i++) 
     68  Y2(i,0) = V2(i,0) + A1(i); 
     69// For each innovation after the first... 
     70for(j=1; j<Nv[1]; j++) { 
     71  // For each time series... 
     72  for(i=0; i<Nv[0]; i++) { 
     73    // The drift plus the current innovation... 
     74    sum = A1(i) + V2(i,j); 
     75    // ...plus the dot product for the VAR(1) term... 
     76    for(k=0; k<Nv[0]; k++) 
     77      sum += B2(i,k) * Y2(k,j-1); 
     78    // ...is the modeled output 
     79    Y2(i,j) = sum; 
    4680  } 
    4781} 

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

    r140 r141  
    2222import os.path 
    2323import scipy as s 
     24from scipy import linalg 
    2425from scipy.stats import distributions 
    2526 
    …… 
    2930from twisted_goodies.pybywire.params import Parameterized 
    3031from twisted_goodies.pybywire.pack import Packer, Unpacker 
     32 
     33import weave 
    3134 
    3235 
    …… 
    106109            lambda _: self.client.registerClasses(*Parameterized.registry)) 
    107110        d.addCallback( 
    108             lambda _: self.client.update('newModel', modelObj)) 
     111            lambda _: self.client.update('newModel', self.modelObj)) 
    109112        d.addCallback(lambda _: setattr(self, 'remoteMode', True)) 
    110113        d.addErrback(self._oops) 
    …… 
    244247 
    245248 
    246 class VAR(Weaver): 
     249class VAR(weave.Weaver): 
    247250    """ 
    248251    Construct me with a sequence of one or more time-series objects containing 
    …… 
    262265        p = len(tsList) 
    263266        tsObj = tsList[0] 
    264         self.y = s.empty(p, len(tsObj)) 
     267        self.y = s.empty((p, len(tsObj))) 
    265268        for k in xrange(len(tsList)): 
    266269            if k > 0: 
    …… 
    269272            self.y[k,:] = tsObj() 
    270273 
    271     def reverse(self, 'a', 'b'): 
    272         return self.c('r', 'y', 'a', 'b', r=s.empty_like(self.y), a=a, b=b) 
     274    def reverse(self, a, b): 
     275        return self.c( 
     276            'v', 'y', 'a', 'b', v=s.empty_like(self.y), a=a, b=b) 
     277 
     278    def forward(self, v, a, b): 
     279        return self.c( 
     280            'y', 'v', 'a', 'b', y=s.empty_like(v), v=v, a=a, b=b) 
    273281 
    274282 
    …… 
    282290    keyAttrs = {'tsList':None} 
    283291    paramNames = ['a', 'b', 'c', 'd', 'e', 'f'] 
    284     support = resource_string(__name__, "%s_support.c" % name) 
    285292 
    286293    def _get_residualizer(self): 
    …… 
    304311        """ 
    305312        """ 
    306  
    307     def  
     313        pass 
    308314 
    309315    def rDistribution(self, N): 

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

    r132 r141  
    1717 
    1818""" 
    19 Unit tests for pyfolio.model 
     19Unit tests for svpmc.model 
    2020""" 
    2121 
    …… 
    2929 
    3030import model 
    31 from data import Data 
    32 from dist.base import Momentum_Distribution 
    33  
     31from tseries import TimeSeries 
    3432import util 
    3533 
    …… 
    8684    def update(self, *args): 
    8785        if args[0] == 'newModel': 
    88             self.model = args[2] 
     86            self.model = args[1] 
    8987        else: 
    9088            raise Exception("Bogus update call!") 
    …… 
    106104class Test_ModelCaller_Basics(util.TestCase): 
    107105    def setUp(self): 
    108         self.models = {1:Mock_Model()} 
     106        self.model = Mock_Model() 
    109107        self.JobClient = model.jobs.JobClient 
    110108        model.jobs.JobClient = Mock_Client 
    111         self.caller = model.ModelCaller(self.models) 
     109        self.caller = model.ModelCaller(self.model) 
    112110 
    113111    def tearDown(self): 
    …… 
    122120        Mock_Model.clearCalls() 
    123121        X = s.array([0.0, 0.1]) 
    124         return self.caller.call(1, 'likelihood', X).addCallback(gotResult) 
     122        return self.caller.call('likelihood', X).addCallback(gotResult) 
    125123 
    126124    def test_setRemoteMode_generic(self): 
    …… 
    131129                'registerClasses', params.Parameterized.registry) 
    132130            self.nextCall( 
    133                 'update', ('newModel', 1, self.models[1])) 
     131                'update', ('newModel', self.model)) 
    134132            self.nextCall( 
    135                 'run', ('callModel', 1, 'likelihood', pack.Packer(X)())) 
     133                'run', ('callModel', 'likelihood', pack.Packer(X)())) 
    136134            for call in Mock_Client.calls: 
    137135                if call[0].endswith('.update') and call[1][0] == 'parallow': 
    …… 
    146144        X = s.array([0.0, 0.1]) 
    147145        d = self.caller.setRemoteMode(None) 
    148         d.addCallback(lambda _: self.caller.call(1, 'likelihood', X)) 
     146        d.addCallback(lambda _: self.caller.call('likelihood', X)) 
    149147        d.addCallback(gotResult) 
    150148        return d 
    …… 
    153151class Test_Prior(util.TestCase): 
    154152    def test_uniformPrior(self): 
    155         p = model.Prior(name='uniform', loc=1, scale=2) 
     153        p = model.Prior(dname='uniform', loc=1, scale=2) 
    156154        X = p.rvs(1000) 
    157155        self.failUnlessEqual(X.shape, (1000,)) 
    …… 
    161159 
    162160 
    163 class Momentum_Normal_Distribution(Momentum_Distribution): 
    164     paramNames = ('sdev',) 
    165     xMax = 6.0 
    166     def H(self, p, t): 
    167         return -(self.sdev**2 * p**2)/2 
    168  
    169  
    170 class Mock_ProjectManager(util.Mock): 
    171     paramNames = ['drift', 'vol'] 
    172     priors = [ 
    173         model.Prior(name='norm', loc=0, scale=1), 
    174         model.Prior(name='uniform', loc=0.5, scale=1.5)] 
    175     expressions = {} 
    176  
    177  
    178 class BaseTC(util.TestCase): 
    179     def makeModel(self): 
    180         self.prices = Data('prices', data=PRICES) 
    181         self.distribution = Momentum_Normal_Distribution() 
    182         self.model = model.Model( 
    183             N_normParams=1, data=self.prices, distribution=self.distribution) 
    184         self.projectManager = Mock_ProjectManager() 
    185         self.projectManager.models = [ 
    186             (self.model, self.projectManager.paramNames)] 
    187         self.mgr = model.ModelManager(self.projectManager) 
    188  
    189  
    190 class Test_Local(BaseTC): 
     161class Test_VAR(util.TestCase): 
    191162    def setUp(self): 
    192         self.makeModel() 
    193  
    194     def test_rPriors(self): 
    195         X = self.mgr.rPriors(1000) 
    196         self.failUnlessAlmostEqual(s.mean(X[:,0]), 0.0,  1) 
    197         self.failUnlessAlmostEqual(s.std(X[:,0]),  1.0,  1) 
    198         self.failUnlessAlmostEqual(s.mean(X[:,1]), 1.25, 1) 
    199         self.failUnless(min(X[:,1]) > 0.5) 
    200         self.failUnless(max(X[:,1]) < 2.0) 
    201  
    202     def test_pPrior(self): 
    203         X = s.array([ 
    204             [-1, 1],    #  
    205             [ 0, 1],    # 1 
    206             [ 0, 0],    #  
    207             [ 1, 1],    # 2 
    208             [ 3, 1],    #  
    209             [ 1, 0],    #  
    210             [ 1, 2],    # 
    211             [ 0, 3],    # 
    212             ]) 
    213         P1 = self.mgr.pPriors(X) 
    214         P2 = s.exp(s.log(P1).sum(0)) 
    215         self.failUnlessAlmostEqual(P2[1], 0) 
    216  
    217     def test_returnsData(self): 
    218         sp = self.fig.add_subplot(111) 
    219         self.model.normParams = [0] 
    220         sp.plot(self.model.normalize(self.model.data())) 
    221  
    222     def test_rDistribution(self): 
    223         N = 10000 
    224         self.model.setParamVector(s.array([0.5, 1.0])) 
    225         RV = self.model.rDistribution(N) 
     163        self.tsList = [ 
     164            TimeSeries('prices', PRICES), TimeSeries('earnings', EARNINGS)] 
     165        self.var = model.VAR(self.tsList) 
     166 
     167    def _plot(self, *vars): 
     168        N = len(vars) 
    226169        fig = self.fig 
    227         sp = fig.add_subplot(111) 
    228         Y, X = s.histogram(RV, bins=50) 
    229         sp.plot(X, Y/(s.diff(X)[0]*N), ls='steps') 
    230         sp.grid(True) 
    231  
    232  
    233 class Test_Remote(BaseTC): 
    234     class CopyableReturner(pb.Root): 
    235         def __init__(self, copyable): 
    236             self.copyable = copyable 
    237         def remote_giveMeCopy(self, null): 
    238             return self.copyable 
    239  
    240     def setUp(self): 
    241         """ 
    242         Create a pb server using L{CopyableReturner} protocol and connect a 
    243         client to it. 
    244         """ 
    245         self.makeModel() 
    246         return self.getReferenceToRoot(self.CopyableReturner(self.model)) 
    247  
    248     def test_remoteVersion(self): 
    249         def got(result): 
    250             self.failIfEqual(id(result), id(self.model)) 
    251             self.failUnless(isinstance(result, model.Model)) 
    252             self.failUnlessEqual(result.N_normParams, self.model.N_normParams) 
    253             self.failUnlessEqual( 
    254                 result.distribution.paramNames, 
    255                 self.model.distribution.paramNames) 
    256             localData = self.model.data() 
    257             remoteData = result.data() 
    258             self.failUnlessElementsEqual(localData, remoteData) 
    259  
    260         params.paregister(*self.model.registry) 
    261         return self.ref.callRemote("giveMeCopy", self.model).addCallback(got) 
    262  
    263  
    264 class Test_ModelAndManager(util.TestCase): 
    265     class ProjectManager(object): 
    266         expressions = {} 
    267         paramNames, priors, models = util.gaussianModelFactory() 
    268      
    269     def setUp(self): 
    270         self.JobClient = model.jobs.JobClient 
    271         model.jobs.JobClient = Mock_Client 
    272         self.projectManager = self.ProjectManager() 
    273  
    274     def tearDown(self): 
    275         model.jobs.JobClient = self.JobClient 
    276  
    277     @defer.deferredGenerator 
    278     def _L(self, *null): 
    279         X = s.linspace(0.1, 4.0, 40) 
    280         L = s.empty_like(X) 
    281         for k, x in enumerate(X): 
    282             wfd = defer.waitForDeferred( 
    283                 self.mgr.likelihoods(s.array([0.0, x]))) 
    284             yield wfd 
    285             L[k] = wfd.getResult() 
    286         I = s.isfinite(L) 
    287         X, L = X[I], L[I] 
    288         sp = self.fig.add_subplot(111) 
    289         sp.plot(X, L, 'o') 
    290         sp.grid(True) 
    291         yield X[L.argmax()] 
    292  
    293     @defer.deferredGenerator 
    294     def test_paramExpressions(self): 
    295         modelObj = self.projectManager.models[0][0] 
    296         for expr in ('sigma+1', 'sigma', '2*sigma'): 
    297             self.projectManager.models = [(modelObj, ['drift', expr])] 
    298             self.mgr = model.ModelManager(self.projectManager) 
    299             wfd = defer.waitForDeferred(self._L()) 
    300             yield wfd; wfd.getResult() 
    301  
    302     def test_pDistribution(self): 
    303         def done(results): 
    304             N_dist = len(results) 
    305             fig = self.fig 
    306             k = 1 
    307             for XY in results.itervalues(): 
    308                 sp = fig.add_subplot(100*N_dist+10+k) 
    309                 sp.plot(XY[0], XY[1], '.') 
    310                 sp.grid(True) 
    311                 k += 1 
    312          
    313         self.mgr = model.ModelManager(self.projectManager) 
    314         d = self.mgr.pDistribution(s.array([0.0, 1.0]), 1000) 
    315         d.addCallback(done) 
    316         return d 
    317  
    318     def test_likelihoods_locally(self): 
    319         self.mgr = model.ModelManager(self.projectManager) 
    320         return self._L().addCallback(self.failUnlessAlmostEqual, 1.0, 1) 
    321      
    322     def test_likelihoods_remotely(self): 
    323         self.mgr = model.ModelManager(self.projectManager) 
    324         d = self.mgr.setRemoteMode(None) 
    325         d.addCallback(self._L) 
    326         d.addCallback(self.failUnlessAlmostEqual, 1.0, 1) 
    327         return d 
    328  
    329          
    330  
     170        for k, var in enumerate(vars): 
     171            sp = fig.add_subplot(100*N+k+11) 
     172            sp.plot(var) 
     173            sp.grid(True) 
     174 
     175    def test_forward_indep(self): 
     176        a = s.array([1.0, 1.0]) 
     177        b = s.array([[0.8, 0.0], [0.0, 0.2]]) 
     178        v = 0.1*s.randn(2, 100) 
     179        y = self.var.forward(v, a, b) 
     180        self._plot(y[0,:], y[1,:]) 
     181 
     182    def test_reverse_indep(self): 
     183        a = s.array([1.0, 1.0]) 
     184        b = s.array([[0.8, 0.0], [0.0, 0.2]]) 
     185        v1 = 0.1*s.randn(2, 100) 
     186        _y = self.var.y 
     187        self.var.y = self.var.forward(v1, a, b) 
     188        v2 = self.var.reverse(a, b) 
     189        self._plot(v1[0,:], v2[0,:], v1[1,:], v2[1,:]) 
     190        self.var.y = _y 
     191 
     192    def test_forward_xcorr(self): 
     193        a = s.array([0.0, 0.0]) 
     194        b = s.array([[0.0, 0.5], [0.0, 0.0]]) 
     195        v = 0.1*s.randn(2, 100) 
     196        y = self.var.forward(v, a, b) 
     197        self._plot(v[0,:], v[1,:], y[0,:]) 
     198 
     199    def test_reverse_xcorr(self): 
     200        a = s.array([1.0, 1.0]) 
     201        b = s.array([[0.01, 0.5], [0.01, 0.01]]) 
     202        v1 = 0.1*s.randn(2, 100) 
     203        _y = self.var.y 
     204        self.var.y = self.var.forward(v1, a, b) 
     205        v2 = self.var.reverse(a, b) 
     206        for k in (0,1): 
     207            self._plot(v1[0,:], v1[1,:], self.var.y[k,:], v2[k,:]) 
     208        self.var.y = _y