Changeset 156

Timestamp:

04/23/08 15:34:26 (8 months ago)

Author:

edsuom

Message:

Working on Model.logLikelihood: VAR process tested

Files:

• projects/AsynCluster/trunk/svpmc/model.c (modified) (5 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (8 diffs)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (5 diffs)

projects/AsynCluster/trunk/svpmc/model.c

@@ -165 +165 @@
-double L = 0;
-
+py::tuple result(2);
+
-
-// For each sample...
@@ -175 +177 @@
-  for(k1=0; k1<Nz[0]; k1++) {
-    // The volatility offset plus the shock...
-2, k1
+1, k0
-    // ...plus the dot product for the VAR(1) term...
-    for(km=0; km<Nz[0]; km++)
-1
+m
-    // ...is the current log-volatility value
-    H2(k1, k0+1) = sum;
-  }
-  
+
-  // Now inverse-scale the innovations for this sample by the square root of
-  // the volatilities in preparation for decorrelating the innovations
@@ -197 +199 @@
-    sum = 0;
-    for(km=0; km<Nz[0]; km++)
-1
+m
-    // ...is the decorrelated multivariate innovation for this sample, given
-    // its log-volatility
@@ -206 +208 @@
-  // multivariate innovation for this sample, conditional on the just-computed
-  // log-volatility h[t].
-    
+
-  // For each time series...
-  for(k1=0; k1<Nz[0]; k1++) {
@@ -223 +225 @@
-    L -= 0.5 * H2(k1, k0+1);
-  }
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+

projects/AsynCluster/trunk/svpmc/model.py

@@ -191 +191 @@
-    _logU_index = -1
-
-n0':500, 'n3
+Mv':100, 'Mz
-    paramNames = params.PARAM_NAMES
-
@@ -242 +242 @@
-        return cv
-
-
-
-
-
+
+
+
+
+
+
+
+
-    def logUniform(self):
-        """
@@ -261 +265 @@
-        return logU[self._logU_index]
-
+    def zProposal(self, z, sigma):
+        """
+        Does a random walk from the supplied 2-D array I{z} of normal random
+        variates, using a symmetric normal proposal distribution with standard
+        deviation I{sigma}.
+        """
+        z = z.copy()
+        self.inline(
+            'z', 'M', 'wiggle', M=self.Mz, wiggle=sigma/s.sqrt(self.Mz))
+        return z
+    
-    def logVolatilities(self, z, v, x, h0):
-        """
@@ -277 +292 @@
-        The method returns a tuple with the total log-likelihood value and the
-        computed log-volatility array I{h}.
+
+        CALLS TO THIS METHOD IS WHERE MOST OF THE CPU TIME IS SPENT!
-        """
-        # Mostly empty log-volatility array
@@ -293 +310 @@
-            self.cache['lv_work'] = y, ri
-        # Run the C code where the heavy lifting is done
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-    def hybridGibbs(self, z, x, sigma):
-        """
@@ -316 +323 @@
-        L_total = 0
-        if 'hg_work' in self.cache:
-v
+
-        else:
-            # Concurrent cross-correlations between volatility shocks
-            rv = s.matrix(
-                linalg.cholesky(self.covarMatrix(self.f), lower=True))
-v
-v
+
+
-        # Initialize volatility shocks from cross-correlation of the IID
-        # variates...
@@ -328 +335 @@
-        # ...and initial log-volatilites, from the offset alone...
-        h0 = self.d.copy()
-        # ...and the log-likelihood of the first sample given the current IID
-        # variates...
-        L = self.logVolatilities(z[:,:N], v[:,:N], x[:,:N], h0)[0]
-        # ...and a set of random-walk proposals for the IID variates, along
-        # with their volatility shocks
@@ -337 +341 @@
-        # Do conditional draws of each sample in turn
-        for k in xrange(self.n):
+            # We need to compare the current and proposed log-volatilities for
+            # this sample, given the previous log-volatility sample
+
+            # Current
+            LV = self.logVolatilities(z[:,k:k+N], v[:,k:k+N], x[:,k:k+N], h0)
+            # Proposal
-            zpk = s.column_stack([zp[:,k], z[:,k+1:k+N]])
-            vpk = s.column_stack([vp[:,k], v[:,k+1:k+N]])
-p, h
+Vp
-            # Metropolis-Hastings accept/reject of the proposal
-
+
+
-                z[k] = zp[k]
-                v[k] = vp[k]
-
-
-
-
+
+
+
-        return L_total
-    

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

@@ -21 +21 @@
-
-import scipy as s
-
+, linalg
-from scipy.signal import signaltools
-
@@ -166 +166 @@
-                "'Uncorrelated' hypothesis not disproved, p=%f" % corrInfo[1])
-    
-    def _check_xcorr(self, dPeak):
-        n = self.model.n
-        for j in xrange(20):
-            h = self.model.draw_h(0.5)
-        # Correlation
-        nd = 6
-        d = s.arange(-nd, nd+1)
-        r = s.array(
-            [s.correlate(h[0,nd+lag:n-nd+lag], h[1,nd:n-nd])[0] for lag in d])
-        lagMetric = r[nd+dPeak] / \
-                    max(s.concatenate([r[:nd+dPeak], r[nd+dPeak:]]))
-        self.failUnless(lagMetric < 10)
-        if VERBOSE:
-            sp = self.fig.add_subplot(111)
-            sp.plot(d, r, 'o-')
-            sp.grid(True)
-    
-
-class Test_VAR(Multivariate_BaseCase):
@@ -263 +246 @@
-
-
-draw_h
+logVolatilities_VAR
-    def setUp(self):
-        self.model = model.Model(tsList=util.TS_LIST)
-
+    def _draw_h(self):
+        self.model.c = s.array([0.0])
+        self.model.g = s.array([0.0])
+        z = s.randn(self.model.p, self.model.n)
+        rv = s.matrix(
+            linalg.cholesky(self.model.covarMatrix(self.model.f), lower=True))
+        v = s.array(rv * z)
+        return self.model.logVolatilities(
+            z, v, s.zeros_like(z), self.model.d.copy())[2]
+    
-    def test_indep_offset(self):
-        n = self.model.n
-        self.model.d = s.array([1.0, 1.0])
-[0.0]
+s.array([0.0])
-        for j in xrange(10):
-            e1, e2 = 0.8*s.rand(2) + 0.1
-            self.model.e = s.array([[e1, 0.0], [0.0, e2]])
-            for k in xrange(20):
-model.draw_h(0.5
+_draw_h(
-            for m, em in enumerate((e1, e2)):
-                ratio = h[m,n/2:].mean() /  ( 1.0/(1-em) )
@@ -286 +279 @@
-        self.model.e = s.array([[e1, 0.0], [0.0, e2]])
-        self.model.f = [0.0]
-
-
+
-        # Uncorrelated
-        self._check_uncorr(h[0,:], h[1,:])
@@ -301 +293 @@
-            sp.plot(h[k,:])
-
+    def _check_xcorr(self, dPeak):
+        n = self.model.n
+        h = self._draw_h()
+        # Correlation
+        nd = 6
+        d = s.arange(-nd, nd+1)
+        r = s.array(
+            [s.correlate(h[0,nd+lag:n-nd+lag], h[1,nd:n-nd])[0] for lag in d])
+        lagMetric = r[nd+dPeak] / \
+                    max(s.concatenate([r[:nd+dPeak], r[nd+dPeak:]]))
+        self.failUnless(lagMetric < 10)
+        if VERBOSE:
+            sp = self.fig.add_subplot(111)
+            sp.plot(d, r, 'o-')
+            sp.grid(True)
+
-    def test_xcorr_var(self):
-        self.model.d = s.array([0.0, 0.0])
-5
+9
-        self.model.f = [0.0]
-        self._check_xcorr(1)