com.meliorbis.numerics.markov.DiscreteStochasticProcessFactory Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of Numerics Show documentation
Show all versions of Numerics Show documentation
A library for working with large multi-dimensional arrays and the functions they represent
package com.meliorbis.numerics.markov;
import static com.meliorbis.numerics.DoubleArrayFactories.createArray;
import static com.meliorbis.numerics.DoubleArrayFactories.createArrayOfSize;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.distribution.NormalDistribution;
import com.meliorbis.numerics.generic.primitives.DoubleArray;
import com.meliorbis.utils.Utils;
/**
* This class has static utility methods for constructing discrete stochastic processes
*/
public class DiscreteStochasticProcessFactory
{
public DiscreteStochasticProcessFactory()
{
}
/**
* Creates an N-state markov transition assuming symmetry around the mean and
* autocorrelation rho using the Rouwenhorst method.
*
* See Kopecky & Suen, 2010.
*
* @param N_ Number of states
* @param rho_ Desired auto-correlation
* @param sigma_epsilon_ The standard deviation of the i.i.d error
*
* @return A pair where the left item contains the N discrete states, and the right item
* the N-state symmetric markov transition matrix with auto-correlation rho
*/
public DiscreteMarkovProcess rouwenhorst(int N_, double rho_, double sigma_epsilon_)
{
double p = (rho_+1d)/2;
// Calculate the upper and lower bounds
double psi = Math.sqrt(((double)(N_-1))*Math.pow(sigma_epsilon_,2)/(1d-Math.pow(rho_,2)));
double[] states = new double[N_];
double step = 2*psi/(N_-1);
states[0] = -psi;
for(int i = 1; i < N_; i++)
{
states[i] = states[i-1]+step;
}
// If there are an odd number of terms
if(N_ % 2 == 1)
{
// then the middle level should be 0, but there will be a small numerical error. remove it.
states[N_/2] = 0;
}
return new DiscreteMarkovProcessImpl(createArray(states),createTransitionMatrix(N_, p, p));
}
/**
* Creates an N-state markov transition assuming symmetry around the mean and
* autocorrelation rho
*
* @param N_ Number of states
* @param rho_ Desired auto-correlation
*
* @return The N-state symmetric markov transition matrix with auto-correlation rho
*/
DoubleArray createSymmetricTransition(int N_, double rho_)
{
double p = (rho_+1d)/2;
return createTransitionMatrix(N_, p, p);
}
/**
* Creates an N-state markov transition matrix using the Rouwenhorst method on parameters
* p, q.
*
* @param N_ The number of states
* @param p_ The parameter p
* @param q_ The parameter q
*
* @return The N-by-N markov transition matrix for the given parameters
*/
DoubleArray createTransitionMatrix(int N_, double p_, double q_)
{
if(p_ <= 0 || p_ >= 1)
{
throw new RuntimeException("p must be in (0,1)");
}
if(q_ <= 0 || q_ >= 1)
{
throw new RuntimeException("q must be in (0,1)");
}
if(N_ <= 0)
{
throw new RuntimeException("The requested number of states must be positive");
}
PolynomialFunction[] pTerms = new PolynomialFunction[N_];
PolynomialFunction[] qTerms = new PolynomialFunction[N_];
PolynomialFunction constant = new PolynomialFunction(new double[]{1d});
// Prepare the individual terms
PolynomialFunction pBase = new PolynomialFunction(new double[]{p_, 1-p_});
PolynomialFunction qBase = new PolynomialFunction(new double[]{1-q_, q_});
pTerms[N_-1] = constant;
qTerms[0] = constant;
// Raise the individidual terms to the appropriate powers
for(int i = 1; i < N_; i++)
{
qTerms[i] = qTerms[i-1].multiply(qBase);
pTerms[N_ - i - 1] = pTerms[N_ - i].multiply(pBase);
}
// Create the transition matrix to fill
DoubleArray transition = createArrayOfSize(N_, N_);
for(int i = 0; i < N_; i++)
{
// Create the polynomial for that row
PolynomialFunction phi_i = pTerms[i].multiply(qTerms[i]);
double[] phi_i_coefs = phi_i.getCoefficients();
// Put the coefficients in the appropriate place
for(int j = 0; j < N_; j++)
{
transition.set(phi_i_coefs[j], i,j);
}
}
return transition;
}
/**
* Creates an N-state markov transition assuming symmetry around the mean and
* autocorrelation rho using the Tauchen method.
*
* See Tauchen, 1986.
*
* @param N_ Number of states
* @param rho_ Desired auto-correlation
* @param sigma_ The standard deviation of the process
* @param omega_ The maximum distance of the process from the mean in either direction, in standard deviations
*
* @return A discrete stochastic process
*/
public DiscreteMarkovProcess tauchen(int N_, double rho_, double sigma_, double omega_)
{
// Calculate the available levels of the process, distributed evenly around the mean (0)
final DoubleArray levels = createArray(Utils.sequence(-omega_ * sigma_, omega_ * sigma_, N_));
// If there are an odd number of terms
if(N_ % 2 == 1)
{
// then the middle level should be 0, but there will be a small numerical error. remove it.
levels.set(0d, N_/2);
}
// The separation between points
final double halfSeparation = omega_ * sigma_ / N_;
// The standard error of the iid shocks
final double sigma_epsilon = Math.sqrt(sigma_*sigma_*(1d-rho_*rho_));
NormalDistribution stdNormal = new NormalDistribution();
// Create a matrix to hold the transition
final DoubleArray transition = createArrayOfSize(N_, N_);
for(int sourceState = 0; sourceState < N_; sourceState++)
{
// This variable keeps track of the last term added because it is what is subtracted from the next term,
// Since they are 2*halfSeparation apart
double lastAdded = 0d;
for(int targetState = 0; targetState < N_ - 1; targetState++)
{
// Add the term for this value
final double newTerm = stdNormal.cumulativeProbability((levels.get(targetState) - rho_ * levels.get(sourceState) + halfSeparation) /
sigma_epsilon);
// Subtract the last added term
double prob = newTerm - lastAdded;
lastAdded = newTerm;
transition.set(prob,sourceState, targetState);
}
// At the higher boundary we have 1 - cumulative distribution to get to a half with below it
transition.set(1d-lastAdded,sourceState, N_ - 1);
}
return new DiscreteMarkovProcessImpl(levels, transition);
}
/**
* Constructs a discrete Markov process with the given levels and transition probabilities
*
* @param levels_ The levels the process can adopt
* @param trans_ The transition probabilities between levels
*
* @return The constructed Markov process
*/
public DiscreteMarkovProcess create(DoubleArray levels_, DoubleArray trans_)
{
return new DiscreteMarkovProcessImpl(levels_, trans_);
}
}