Changeset 158

Timestamp:

04/25/08 00:24:02 (7 months ago)

Author:

edsuom

Message:

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

Files:

• projects/AsynCluster/trunk/doc/gevolver (deleted)
• projects/AsynCluster/trunk/doc/svpmc (added)
• projects/AsynCluster/trunk/doc/svpmc/L_for_log-volatility_sim.pdf (added)
• projects/AsynCluster/trunk/svpmc/model.c (modified) (2 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (8 diffs)
• projects/AsynCluster/trunk/svpmc/params.py (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (2 diffs)

projects/AsynCluster/trunk/svpmc/model.c

@@ -95 +95 @@
-double norm, normp, L, Lp, dL;
-
+// Seed the PRNG on the first run only
+static int seeded = 0;
+if (!seeded) {
+  int i;
+  int j = 0;
+  unsigned short seed[3];
+  for(i=0;i<3;i++) {
+    do {
+      j++;
+    } while (j < (i+1)*10000);
+    seed[i] = clock()*clock() % 65537;
+  }
+  seed48(seed);
+  seeded = 1;
+}
-
-// for each multivariate sample...
@@ -113 +128 @@
-    for(k2=0; k2<M; k2++) {
-      // Generate a uniform, symmetric, random walk proposal
-
+ 
-      // Do the ratio in log space
-      Lp = -0.5 * pow(normp, 2);

projects/AsynCluster/trunk/svpmc/model.py

@@ -128 +128 @@
-                paramContainer.pLogLikelihood = -s.inf
-            up = pack.Unpacker(result)
-pLogLikelihood
-latent
+L
+z
-
-        def gotLikelihood(result):
-pLogLikelihood
-latent
+L
+z
-        
-        if (remote or self.remoteMode) and not local:
@@ -150 +150 @@
-    Construct me with a 2-D array of time series values, one time series per
-    row.
-
+
-    Then you can call me with a C{p x 1} drift vector I{a} and a C{p x p}
-    VAR(1) cross-correlation matrix I{b}, to get a C{p x n} vector of residuals
@@ -191 +191 @@
-    _logU_index = -1
-
-100, 'Mz':10
+50, 'Mz':5
-    paramNames = params.PARAM_NAMES
-
@@ -250 +250 @@
-        return cv
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    
-    def logUniform(self):
@@ -321 +329 @@
-        return L0, L, h[:,1:]
-
-N=100
+tol=1E-4
-        """
-        Does one iteration of a hybrid Gibbs sampler for the log-volatilities,
@@ -328 +336 @@
-        Returns the likelihood of the innovations in I{x} given the simulated
-        log-volatilities. Updates the IID normal variates in place.
+
+        The method will account for the effect of log-volatility proposals on
+        downstream samples. You can specify the maximum fractional effect to
+        account for with I{tol}. Smaller tolerances require more post-sample
+        log-volatility recomputations and hence more CPU time.
-        """
-        L_total = 0
-        acceptances = 0
-        rv = self.rv
+        N = self.decaySamples(tol)
-        # Initialize volatility shocks from cross-correlation of the IID
-        # variates...
@@ -369 +383 @@
-        Call this method with a parameter object that defines:
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
-
-        A new 2-D array of log-volatilities will be drawn from a hybrid Gibbs
@@ -400 +416 @@
-        self.setParamSequence(paramContainer.paramSequence())
-        # ...including the latent parameter of IID normal variates
-latent
+z
-        # Innovations as residuals of the reversed VAR(1)
-        x = self.var.reverse(self.a, self.b)

projects/AsynCluster/trunk/svpmc/params.py

@@ -239 +239 @@
-    """
-    """
-latent':[], 'pLogLikelihood':None, 'pLogProposal
+z':[], 'L
-    paramNames = PARAM_NAMES
-
-    def __add__(self, other):
-latent=self.latent
+z=self.z
-        newVersion.setParamSequence([
-            getattr(self, name) + getattr(other, name)

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

@@ -288 +288 @@
-        if VERBOSE:
-            self.fig.add_subplot(111).plot(z1, z2)
+
+    def test_decaySamples(self):
+        x = s.zeros((5, 1000))
+        x[2,0] = 1.0 # The impulse
+        var = model.VAR(x)
+        a = s.zeros(5)
+        modelObj = model.Model(tsList=util.TS_LIST)
+        for j in xrange(20):
+            # Keep all coefficients small enough to ensure stability
+            b = modelObj.e = 0.3*s.rand(5,5)
+            y = var.forward(x, a, b)
+            for tol in s.logspace(-5, -1, 10):
+                N = modelObj.decaySamples(tol)
+                print tol, N
+                failed = abs(y[:,N].max()) > tol
+                if (N > 100 or failed) and VERBOSE:
+                    N_plot = min([3*N, 30])
+                    sp = self.fig.add_subplot(111)
+                    sp.semilogy(y[:,:N_plot].transpose())
+                    sp.axvline(N)
+                    sp.axhline(tol)
+                    self.failIf(
+                        failed,
+                        "Estimate of %d samples " % N +\
+                        "to decay within %f tolerance was insufficient" % tol)
-        
-
@@ -507 +532 @@
-
-
-class Test_Model_hybridGibbs(Model_BC_Mixin, Multivariate_BaseCase):
-    N = 10000
-    corr = 0.0
-    params = {
-        'a': [0],
-        'b': [[0]],
-        'c': [],
-        'd': [0.01],
-        'e': [[0.0]],
-        'f': [],
-        'g': [corr],
-        }
-
-
-class Test_Model_likelihood(Model_BC_Mixin, Multivariate_BaseCase):
-5
+30
-    corr = 0.0
-    params = {