Changeset 157

Timestamp:

04/24/08 20:15:13 (9 months ago)

Author:

edsuom

Message:

Got log-volatility Metropolis-within-Gibbs simulation working!!!

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) (8 diffs)

projects/AsynCluster/trunk/svpmc/model.c

@@ -93 +93 @@
-
-int k0, k1, k2;
-
+, dL
-
-
@@ -103 +103 @@
-  // which we build a multivariate log-volatility proposal for this
-  // multivariate sample.
-
+
-  // For each time series...
-  for(k1=0; k1<Nz[0]; k1++) {
@@ -113 +113 @@
-    for(k2=0; k2<M; k2++) {
-      // Generate a uniform, symmetric, random walk proposal
-
+(double)
-      // Do the ratio in log space
-      Lp = -0.5 * pow(normp, 2);
@@ -155 +155 @@
-// The multivariate scratchpad y is used as follows, by column:
-// ----------------------------------------------------------------------------
-2*(1-sigma**2) for bivariate normal probability density function
-0.5*log(pi*2*sigma**2) offset for bivariate normal pdf
+Bivariate norm PDF: sqrt(1-sigma**2)
+Bivariate norm PDF, log-scaling: log(sqrt(2*pi)*sqrt(1-sigma**2))
-// 2    Innovation values after inverse-scaling by log-volatility
-// 3    Inverse-scaled innovation values after decorrelation
@@ -218 +218 @@
-    // The log-likelihood for this sample is simply the argument of the
-    // exponential of the normal pdf, with the reciprocal scaling of the
-
-
+
+
+
-    // Since the innovation has been inverse-scaled by the square root of the
-    // proposed volatility, the probability density needs to be scaled as

projects/AsynCluster/trunk/svpmc/model.py

@@ -191 +191 @@
-    _logU_index = -1
-
-
+0
-    paramNames = params.PARAM_NAMES
-
@@ -223 +223 @@
-    var = property(_get_var)
-
+    def _get_rv(self):
+        if 'rv' not in self.cache:
+            # Concurrent cross-correlations between volatility shocks
+            self.cache['rv'] = s.matrix(
+                linalg.cholesky(self.covarMatrix(self.f), lower=True))
+        return self.cache['rv']
+    rv = property(_get_rv)
+
-
-    #--- Support --------------------------------------------------------------
@@ -273 +281 @@
-        z = z.copy()
-        self.inline(
-/s.sqrt(self.Mz)
+
-        return z
-    
@@ -296 +304 @@
-        """
-        # Mostly empty log-volatility array
-empty
+zeros
-        # Working arrays with intialized values
-        if 'lv_work' in self.cache:
@@ -303 +311 @@
-            # Scratchpad with some constants
-            y = s.empty((self.p, 4))
-2 * 
-0.5 * s.log(s.pi*
+s.sqrt
+s.log(s.sqrt(2*s.pi) * 
-            # Inverse of concurrent cross-correlations between innovation shocks
-            ri = linalg.inv(
@@ -313 +321 @@
-        return L0, L, h[:,1:]
-
-
+, N=100
-        """
-        Does one iteration of a hybrid Gibbs sampler for the log-volatilities,
@@ -322 +330 @@
-        """
-        L_total = 0
-
-
-
-
-
-
-
-
+
+
-        # Initialize volatility shocks from cross-correlation of the IID
-        # variates...
-rv * z
+s.array(rv * z)
-        # ...and initial log-volatilites, from the offset alone...
-.copy()
+
-        # ...and a set of random-walk proposals for the IID variates, along
-        # with their volatility shocks
-        zp = self.zProposal(z, sigma)
-rv * zp;
+s.array(rv * zp)
-        # Do conditional draws of each sample in turn
-        for k in xrange(self.n):
-
-
-
-
+
+
-            LV = self.logVolatilities(z[:,k:k+N], v[:,k:k+N], x[:,k:k+N], h0)
-P
+...to a log-volatility p
-            zpk = s.column_stack([zp[:,k], z[:,k+1:k+N]])
-            vpk = s.column_stack([vp[:,k], v[:,k+1:k+N]])
-            LVp = self.logVolatilities(zpk, vpk, x[:,k:k+N], h0)
-            # Metropolis-Hastings accept/reject of the proposal
-
+
+
+
-                LV = LVp
-k] = zp[
-k] = vp[
+:,k] = zp[:,
+:,k] = vp[:,
-            # Update stuff for the next sample
-
+:,
-            L_total += LV[0]
-
+, float(acceptances)/self.n
-    
-
@@ -401 +403 @@
-        # Innovations as residuals of the reversed VAR(1)
-        x = self.var.reverse(self.a, self.b)
-
-
-
+
+
+
+

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

@@ -246 +246 @@
-
-
+class Test_Model_misc(util.TestCase):
+    def test_logUniform(self):
+        N = 35000
+        import timeit
+        tObserved = timeit.Timer(
+            "M.logUniform()",
+            "import model\nM = model.Model()").timeit(N)
+        tBaseline = timeit.Timer(
+            "s.log(s.rand())",
+            "import scipy as s").timeit(N)
+        percent = 100 * tObserved/tBaseline
+        if VERBOSE:
+            print "The efficient method took %2d%% as long as the stupid way" \
+                  % percent
+        self.failUnless(percent < 40)
+
+    def _check_normality(self, x):
+        p = stats.normaltest(x)[1]
+        if VERBOSE:
+            fig = self.fig
+            sp = fig.add_subplot(211)
+            sp.plot(x)
+            sp = fig.add_subplot(212)
+            sp.hist(x, bins=20)
+        self.failUnless(
+            p > 0.05, "Deviation from normality with significance p=%f" % p)
+
+    def test_zProposal(self):
+        N = 1000
+        z = s.zeros((2,3))
+        modelObj = model.Model()
+        z1 = s.empty(N)
+        z2 = s.empty(N)
+        for k in xrange(N):
+            z = modelObj.zProposal(z, 0.5)
+            z1[k] = z[0,0]
+            z2[k] = z[1,1]
+        self._check_normality(z1[0.7*N:])
+        self._check_normality(z2[0.7*N:])
+        self.failUnless(stats.pearsonr(z1, z2)[1] > 0.05)
+        if VERBOSE:
+            self.fig.add_subplot(111).plot(z1, z2)
+        
+
-class Test_Model_logVolatilities_VAR(Multivariate_BaseCase):
-    def setUp(self):
@@ -252 +296 @@
-    def _draw_h(self):
-        self.model.c = s.array([0.0])
-
+, 0.0
-        z = s.randn(self.model.p, self.model.n)
-        rv = s.matrix(
@@ -259 +303 @@
-        return self.model.logVolatilities(
-            z, v, s.zeros_like(z), self.model.d.copy())[2]
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    def test_indep_offset(self):
-        n = self.model.n
@@ -272 +338 @@
-                ratio = h[m,n/2:].mean() /  ( 1.0/(1-em) )
-                self.failUnlessAlmostEqual(ratio, 1.0, 0)
-
+
-    def test_indep_shocks(self):
-        n = self.model.n
@@ -322 +388 @@
-
-
-
+logVolatilities_
-    def setUp(self):
-        self.model = model.Model()
-
-
-
-
+
+
-            # Tweak the model object to make it testable for this p value
-            self.model.tsList = [util.TS_LIST[0]] * p
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-                # Check for non-correlation
-                for j in xrange(p):
@@ -347 +421 @@
-                        if j == k:
-                            continue
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-    def test_complex(self):
-
-
+
+
+
-            # Stuff specific to this number of time series
-            N_correlations = sum([p-k-1 for k in xrange(p)])
-            cMax = 0.90 / s.sqrt(p/1.5)
-            # Tweak the model object to make it testable for this p value
-            self.model.tsList = [util.TS_LIST[0]] * p
-            for i in xrange(10):
-                # Randomly generate a set of cross-correlations and set the
-                # model's innovation shock correlations parameter
-self.model.
+
-                # Generate some multivariate normal test data having the
-                # generated correlations
-
-
-
+
+
+
+
-                x = random.multivariate_normal(s.zeros(p), r, n).transpose()
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    N = 10000
-    corr = 0.0
@@ -389 +519 @@
-        'g': [corr],
-        }
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-    class Mock_ParameterContainer(object):
@@ -403 +545 @@
-        return pc
-    
-
-
+
+
+
+
-        # Simulate a single time series of data per the parameters
-        self.iv = random.multivariate_normal(
-
-
-
-
-
-
-
+
+
+
+
+
+
-        # Instantiate a model to be tested and set its parameters
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        if VERBOSE:
-            fig = self.fig
-iv[1,:], self.x, p, z[0
+x[0,:], self.h, h[0,:
-            for k, vector in enumerate(vectors):
-
-
-
+
+
+
+
+
-        
+    def test_hybridGibbs_highCorrelation(self):
+        self._setupModel(e=[[0.95]])
+        # Sigma=0.5 seems to give us the supposedly ideal ~44% acceptance rate
+        self._runHG(50, 0.5)
+        
+    def test_hybridGibbs_lowCorrelation(self):
+        self._setupModel(e=[[0.4]])
+        self._runHG(50, 0.5)
+
+    def test_hybridGibbs_noCorrelation(self):
+        self._setupModel(e=[[0.0]])
+        self._runHG(50, 0.5)