import numpy
import scipy

def fast_sample(w, M, u):
    """Generates a sample of size M of integers with weights given by w.
    
    Parameters
    ----------
    w : array_like
        A vector of real weights corresponding to each of the integers ``0`` to
        ``(w.size - 1)``.  These weights do not need to be normalized.
    M : integer
        The requested sample size.
    u : real
        A Uniform(0,1) draw, used as an offset to the CDF.
    
    Returns
    -------
    sample : array_like
        An array containing the resulting sample.
    """

    # Normalize the weights so that they sum to 1.
    w = w / sum(w)

    # Initialize storage for the cumulative weight CDF (including the ``u``
    # offset) and the resulting sample.
    c = numpy.empty((M))
    sample = numpy.empty((M), dtype=numpy.int)

    # Draws for the first integer (0).  a and b are the beginning and ending
    # offsets for each integer.
    c[0] = M * w[0] + u
    a = 0
    b = numpy.math.floor(c[0])
    if b > a:
        sample[a:b] = 0

    for j in numpy.arange(1, w.size):
        c[j] = c[j-1] + M * w[j]
        a = b
        b = numpy.math.floor(c[j])
        if b > a:
            sample[a:b] = j
    return sample

def _test():
    seed = 274
    numpy.random.seed(274)

    N = 10
    M = 1000
    u = numpy.random.rand()
    w = numpy.random.rand(N)
    sample = fast_sample(w, M, u)
    num, bins = scipy.histogram(sample, numpy.arange(N))
    u = 0.2
    print 'u = ', u
    print 'w = ', w / sum(w)
    print 'frequency = ', num / float(M)
    print 'w - sample frequency = ', w / sum(w) - num / float(M)

    u = 0.98
    sample = fast_sample(w, M, u)
    num, bins = scipy.histogram(sample, numpy.arange(N))
    print 'u = ', u
    print 'w = ', w / sum(w)
    print 'frequency = ', num / float(M)
    print 'w - sample frequency = ', w / sum(w) - num / float(M)

if __name__ == "__main__":
    _test()
Content-Type: text/plain; charset=UTF-8
Content-Length: 1961
Content-Disposition: inline; filename="ou_gdp.py"
Last-Modified: Wed, 19 Now 2008 12:02:31 GMT
Expires: Wed, 19 Now 2008 12:07:31 GMT

# ou-gdp.py
#
# Estimates an Ornstein-Uhlenbeck process using data on U.S. GDP.
#
# Jason Blevins <jrblevin@sdf.lonestar.org>
# Durham, March 24, 2008 09:02 EDT

from numpy import *
from scipy import stats, optimize
from pylab import *


def prt(x):
    print "theta = %f\t%f" % (x[0], x[1])

## Log-likelihood function
def log_likelihood(theta, t, x):
    eta = theta[0]
    sigma = exp(theta[1])                         # sigma >= 0

    result = 0.0

    # Calculate the likelihood, skipping the first observation.
    for i in arange(1, size(t)):

        # Time interval (here, constant)
        dt = t[i] - t[i-1]

        # Standard deviation of x(i)
        sd = sigma * sqrt((1 - exp(-2*eta*dt)) / (2 * eta))

        # Mean of x(i) (x_bar is zero since we de-meaned already)
        mean = x[i-1] * exp(-eta * dt)

        # Accumulate the result
        result += log(stats.norm.pdf(x[i], loc=mean, scale=sd))

    return result / size(t)

# Paths
data_dir = '/home/jrblevin/projects/macro-data'
gdp_file = data_dir + '/us-gdp.dat'
deflator_file = data_dir + '/us-gdp-deflator.dat'

# Load the data
yr, qtr, gdp = loadtxt(gdp_file, unpack=True)
deflator = loadtxt(deflator_file, usecols=[2])
log_gdp = log(gdp)

# Linear regression
t = arange(size(log_gdp))
(alpha, beta) = polyfit(t, log_gdp, 1)
fit = polyval([alpha, beta], t)
resid = log_gdp - alpha * t - beta

# Plot the data and the fitted line
# plot(log_gdp, label="Log Real GDP")
# plot(fit, label="Fitted line")
# title('Log Real GDP Regression')
# show()

# Starting values (a nonnegative transformation is applied to sigma)
eta = 0.1
sigma = 0.1
theta = (eta, log(sigma))

# Define a pure function which returns the negative log-likelihood
# at the given parameter values.
f = lambda theta: -log_likelihood(theta, t, resid)

theta = optimize.fmin(f, theta)

print "eta = %f" % theta[0]
print "sigma = %f" % exp(theta[1])

# Plot the residuals
#plot(resid)
#title('GDP Residuals')
#show()
Content-Type: text/plain; charset=UTF-8
Content-Length: 1168
Content-Disposition: inline; filename="ou_sim.py"
Last-Modified: Wed, 19 Now 2008 12:02:31 GMT
Expires: Wed, 19 Now 2008 12:07:31 GMT

# ou_sim.py
#
# Simulates an Ornstein-Uhlenbeck process.
#
# Jason Blevins <jrblevin@sdf.lonestar.org>
# Durham, March 24, 2008 09:16 EDT

import numpy
from numpy import exp, sqrt

def ou_sim(T, sigma=1.0, eta=1.0, mu=0.0):
    """Ornstein-Uhlenbeck Simulator

    Parameters
    ----------

    T -- number of periods to simulate

    sigma -- standard deviation of the Brownian Motion
    eta -- speed of mean reversion
    mu -- mean

    Returns
    -------

    x -- an array containing the realization of the process
    """
    x = numpy.ones((T)) * mu                  # Initialize storage for x

    # Store a series of independent Normal(0,1) draws
    omega = numpy.random.normal(size=T)

    for t in numpy.arange(1,T):
        x[t] = x[t-1] * exp(-eta) + mu * (1 - exp(-eta)) \
            + sigma * sqrt((1 - exp(-2*eta))/(2*eta)) * omega[t]

    return x

if __name__ == '__main__':
    import pylab

    sigma = 0.01
    eta = 0.03
    mu = 0.0

    x = ou_sim(250, sigma, eta, mu)

    pylab.plot(x, label = ("eta = %0.2f, sigma=%0.2f" % (eta, sigma)))
    pylab.title('Ornstein-Uhlenbeck process')
    pylab.legend(loc='upper left')
    pylab.show()