com.meliorbis.economics.individual.egm.EGMSolverBase Maven / Gradle / Ivy
package com.meliorbis.economics.individual.egm;
import static com.meliorbis.numerics.DoubleArrayFactories.createArrayOfSize;
import static com.meliorbis.numerics.generic.primitives.impl.Interpolation.interpolateFunction;
import static com.meliorbis.numerics.generic.primitives.impl.Interpolation.interpolateFunctionAcross;
import java.util.LinkedList;
import java.util.List;
import org.apache.commons.lang.ArrayUtils;
import com.meliorbis.economics.individual.IndividualProblemSolver;
import com.meliorbis.economics.infrastructure.SolverBase;
import com.meliorbis.economics.model.Model;
import com.meliorbis.economics.model.ModelConfig;
import com.meliorbis.economics.model.ModelException;
import com.meliorbis.economics.model.State;
import com.meliorbis.numerics.generic.primitives.DoubleArray;
import com.meliorbis.numerics.generic.primitives.DoubleIndexedReduction;
import com.meliorbis.numerics.generic.primitives.DoubleSubspaceSplitIterator;
import com.meliorbis.numerics.generic.primitives.impl.DoubleArrayFunctions;
import com.meliorbis.numerics.generic.primitives.impl.Interpolation.Params;
import com.meliorbis.numerics.threading.ComputableRecursiveAction;
import com.meliorbis.utils.Timer;
import com.meliorbis.utils.Timer.Stoppable;
import com.meliorbis.utils.Utils;
/**
* The base class for solvers that use the method of endogenous gridpoints
*
* @author Tobias Grasl
*
* @param The config class
* @param The state class
* @param The model class
*/
public abstract class EGMSolverBase,
M extends Model> extends SolverBase
implements IndividualProblemSolver
{
protected DoubleArray _normalisedSOPStates;
protected DoubleArray _eopStates;
protected Params _interpParams;
protected int[] _futureConditionalContribsSize;
protected int _firstAggStateDim;
protected int _firstIndividualStateDim;
protected int[] _arrangeIndexs;
protected int[] _acrossDims;
private int _nIndExpn = 1;
/**
* Constructs the solver for the given model and configuration
*
* @param model_ The model this solver is for
* @param config_ The configuration to apply
*/
public EGMSolverBase(M model_, C config_)
{
super(model_, config_);
}
/**
* This default implementation does nothing, override if other behaviour is necessary
*
* @param state_ The calculation state
*/
@Override
public void initialise(S state_)
{
state_.setIndividualError(Double.MAX_VALUE);
}
/**
* Creates an array appropriate for holding the conditional future contributions
*
* @return The created array
*/
protected DoubleArray createFutureConditionalContributionsGrid()
{
return createArrayOfSize(_futureConditionalContribsSize);
}
/**
* Sets the end-of-period individual state values which correspond to the start-of-period state values on the
* endogenous grid. This allows interpolation from e-o-p to s-o-p states.
*
* The array must be of individual variable grid size as returned by {@link Model#createIndividualVariableGrid()}.
*
* @param eopStates_ The individual state values
*/
final protected void setEndOfPeriodStates(DoubleArray eopStates_)
{
_eopStates = eopStates_;
}
/**
* Sets the start-of-period individual state values which form the on-grid 'endogenous states' used during the calculation
* of conditional future expectations in {@link EGMSolverBase#calculateFutureConditionalExpectations(int[], int[], State)}.
*
* The array must be of individual transition grid size as returned by {@link Model#createIndividualTransitionGrid()}. The states must be
* normalised by the normalising aggregate shock and any growth parameters.
*
* @param sopStates_ The individual state values
*/
final protected void setStartOfPeriodStates(DoubleArray sopStates_)
{
_normalisedSOPStates = sopStates_;
}
/**
* Each iteration consists of the following steps:
*
* - Interpolate expected future individual state choices (i.e. x''). Performed in {@link #interpExpectedFutureTransition(State)}.
* - Calculate future conditional expectations of individuals in the euler equation. Performed in {@link #calculateFutureConditionalContributions(State)}.
* - Determine unconditional expectations from the conditional ones. Performed in {@link #individualUnconditionalExpectation(DoubleArray, State)}.
* - Calculate the implied start-of-period state in the current period from the expectations. Performed in {@link #calculateImpliedStartOfPeriodState(DoubleArray, State)}.
* - Interpolate those implied states back to the grid. Performed in {@link #interpIndividualTransition(DoubleArray, State)}.
*
* The second (indirectly) and fourth steps must be implemented by model specific code, the others are performed automatically but
* can be overridden.
*
* @param state_ The current state of the calculation, which will be updated by this method.
*/
@Override
public final void performIteration(S state_) throws ModelException
{
beginIteration(state_);
// First, determine the expected future choices from the current ones
// calculates so far
interpExpectedFutureTransition(state_);
// First, calculate the individual future conditional expectations for
// any possible individual current end-of-period state
// for all possible current aggregate endogenous vars and all possible
// combinations of current and future exogenous states
DoubleArray futureConditionalContributions = calculateFutureConditionalContributions(state_);
Timer timer = new Timer();
Stoppable stoppable = timer.start("expectation");
// Collapse the individual conditional expectations to expectations
// given the current state
DoubleArray expectedFutureCondition = individualUnconditionalExpectation(futureConditionalContributions, state_);
stoppable.stop();
validateExpectedFutureCondition(expectedFutureCondition, state_);
stoppable = timer.start("calcPeriodStart");
// Use the expectations of future variables given current end-of-period
// states to calculate the implied current start-of-period
// states
final DoubleArray impliedStartOfPeriodState = calculateImpliedStartOfPeriodState(expectedFutureCondition, state_);
stoppable.stop();
validateStartOfPeriodState(impliedStartOfPeriodState, state_);
stoppable = timer.start("interp");
DoubleArray gridValues = interpIndividualTransition(impliedStartOfPeriodState, state_);
stoppable.stop();
DoubleArray oldPolicy = state_.getIndividualPolicy();
// Store the freshly calculated individual transitions
state_.setIndividualPolicy(gridValues);
updateError(oldPolicy, gridValues, state_);
afterPolicyUpdate(oldPolicy, gridValues, state_);
}
/**
* Calculates the error from old and new policy.
* Default calculates each period and compares every value. The error is
* stored back in the state.
*
* @param oldPolicy_ The old policy
* @param newPolicy_ The newly determined policy
* @param state_ The state of the calculation
*/
protected void updateError(DoubleArray oldPolicy_, DoubleArray newPolicy_, S state_)
{
state_.setIndividualError(DoubleArrayFunctions.maximumRelativeDifferenceSpecial(oldPolicy_, newPolicy_));
}
/**
* A callback method called at the beginning of each iteration. Override to perform necessary
* per-iteration pre-processing.
*
* @param state_ The current state of the calculation
*/
protected void beginIteration(S state_)
{
}
/**
* A callback method called when the individual policy function has been recalculated using this
* solver. Override to perform post-processing.
*
* @param oldPolicy_ The old individual state policy.
* @param newPolicy_ The new individual state policy.
* @param state_ The current state of the calculation.
*/
protected void afterPolicyUpdate(DoubleArray oldPolicy_, DoubleArray newPolicy_, S state_)
{
}
/**
* This method interpolates the individual state transition grid along the
* aggregate state dimensions, since the future expected aggregate state is
* not on the grid. Sets the result on {@code state_}.
*
* @param state_ The current state of calculation
*/
protected void interpExpectedFutureTransition(final S state_)
{
state_.setExpectedIndividualTransition(_model.conditionalExpectation(state_.getIndividualPolicy(), state_));
}
/**
* Given the implied start of period states, interpolates to discover the individual transition function
*
* @param impliedInitialStates_ The array of states
* @param state_ The current state of calculation
*
* @return The array representing the function values of the transition from
* the standard grid
*/
protected DoubleArray interpIndividualTransition(final DoubleArray impliedInitialStates_, final S state_)
{
return interpolateIndividualChoiceToGrid(impliedInitialStates_, _eopStates);
}
/**
* Given a mapping from implied initial individual states to a variable, calculates the values for that mapping at
* on-grid start of period states.
*
* @param impliedInitialStates_ The initial states, which can be off grid
* @param choiceVariable_ The function values
*
* @return The function values at on grid points.
*/
protected DoubleArray interpolateIndividualChoiceToGrid(final DoubleArray impliedInitialStates_, DoubleArray choiceVariable_)
{
if (_config.getAggregateNormalisingExogenousStates().size() > 0)
{
// -1 because the last dimension is variables dimension
int[] normDimensions = Utils.sequence(_normalisedSOPStates.numberOfDimensions()
- _config.getAggregateNormalisingExogenousStates().size() - 1, _normalisedSOPStates.numberOfDimensions()-1);
return interpolateFunctionAcross(impliedInitialStates_, choiceVariable_, 0, _normalisedSOPStates,
_interpParams , normDimensions);
} else
{
return interpolateFunction(impliedInitialStates_, choiceVariable_, 0, _normalisedSOPStates,
_interpParams);
}
}
/**
* Given expectations of a variable conditional on each possible joint
* transition of aggregate and individual exogenous states, calculates the unconditional expectations given only the
* current exogenous states..
*
* @param conditionalExpectation_ The conditional expectations from which to determine the unconditional one
* @param state_ The current calculation state
*
* @return The unconditional expectations calculated
*/
protected DoubleArray individualUnconditionalExpectation(DoubleArray conditionalExpectation_, final S state_)
{
DoubleArray transitionProbabilities = _config.getExogenousStateTransition();
DoubleArray results = _model.createIndividualVariableGrid(getNumberOfIndividualExpectations());
Timer timer = new Timer();
// Multiply by the transition probs and sum across the target dims
Stoppable stoppable = timer.start("call unconditional");
final int nIndShocks = _config.getIndividualExogenousStateCount();
final int nIndStates = _config.getIndividualEndogenousStateCount();
// Iterate across exogenous individual states - since these are not in the conditonal
// contributions array because the current ind state does not matter when evaluating potential
// next period outcomes
DoubleSubspaceSplitIterator currentExoIter = transitionProbabilities.iteratorAcross(Utils.sequence(0, nIndShocks));
final int[] multAcrossDims = Utils.sequence(0, transitionProbabilities.numberOfDimensions()-nIndShocks);
final int[] sumAcrossDims = Utils.sequence(1, transitionProbabilities.numberOfDimensions()-nIndShocks);
List actions = new LinkedList<>();
while(currentExoIter.hasNext()) {
currentExoIter.nextDouble();
final int[] currentExoIndex = currentExoIter.getIndex();
actions.add(()->{
// The transition probabilities conditional on the prior individual state
DoubleArray probsForState = transitionProbabilities.at(currentExoIndex);
// Multiply each value with its probability
DoubleArray weighted = conditionalExpectation_.across(multAcrossDims).multiply(probsForState);
// take a sum across the future outcomes
DoubleArray uncondExpForState = (DoubleArray) weighted.across(sumAcrossDims).sum();
final int[] fillAt = Utils.repeatArray(-1, results.numberOfDimensions());
// Copy the individual shocks this calc is for into the index indicating where results are to be stored
System.arraycopy(currentExoIndex, 0, fillAt, nIndStates, nIndShocks);
// copy the results over, arranging the dimensions appropriately
if( _acrossDims != null)
{
results.at(fillAt).fillDimensions(uncondExpForState.arrangeDimensions(_arrangeIndexs),_acrossDims);
}
else
{
results.at(fillAt).fill(uncondExpForState.arrangeDimensions(_arrangeIndexs));
}
});
}
getNumerics().getExecutor().executeAndWait(actions);
stoppable.stop();
return results;
}
/**
* Calculates the expectational part of the Euler equation(s).
*
* This method performs generic, technical parts of determining what must be calculate and calling
* those operations in a multi-threaded way. Actual model-specific calculations are delegated to
* {@link #calculateFutureConditionalExpectations}.
*
* @param state_ The current calculation state.
*
* @return The calculated expectations
*
* @throws ModelException If an error occurs calculating the expectations.
*/
protected final DoubleArray calculateFutureConditionalContributions(S state_) throws ModelException
{
/*
* Create the array to hold the results; multi-threaded filling is not a
* problem
*/
final DoubleArray futureConditionalExpectations = createFutureConditionalContributionsGrid();
futureConditionalExpectations.across(_firstIndividualStateDim, _futureConditionalContribsSize.length - 1).reduce(new DoubleIndexedReduction()
{
@Override
public double perform(DoubleSubspaceSplitIterator iter_) throws ModelException
{
// Get the full index, including the placeHolders!
int[] index = iter_.getFullIndex();
// Get the stochastic shock transition index
int[] transIndex = ArrayUtils.subarray(index, 0, _firstAggStateDim);
// Construct a full transition array which includes the initial index
transIndex = ArrayUtils.addAll(new int[]{0}, transIndex);
// Get the aggregate state index
int[] aggStateIndex = ArrayUtils.subarray(index, _firstAggStateDim, _firstIndividualStateDim);
// Get the probability of that shock transition
//double probability = _config.getExogenousStateTransition().get(transIndex);
// Don't select the dimension across variables because that
// one is also filled
index[index.length - 1] = -1;
index[_firstIndividualStateDim] = -1;
/*
* Ignore 0 probability events
* CAN'T DO THIS BECAUSE OF OPTIMISATION WHICH IGNORES PRIOR IND STATE
*/
//if (probability != 0d)
{
// calculate and store the values
final DoubleArray transitionResults = calculateFutureConditionalExpectations(transIndex, aggStateIndex, state_);
futureConditionalExpectations.at(index).fill(transitionResults);
}
return Double.NaN;
}
});
afterFutureConditionalExpectations(futureConditionalExpectations, state_);
return futureConditionalExpectations;
}
/**
* Callback method for further processing after all the conditional future expectations have been calculated in each
* iteration.
*
* The default method does nothing
*
* @param futureConditionalExpectations_ The freshly-baked future conditional expectations
* @param state_ The state of the ongoing calculation
*/
protected void afterFutureConditionalExpectations(DoubleArray futureConditionalExpectations_, S state_)
{
// Default does nothing
}
/**
* @return The expected dimensionality of the function returned by {@link #calculateFutureConditionalExpectations(int[], int[], State)}. Defaults
* to 1.
*/
final protected int getNumberOfIndividualExpectations()
{
return _nIndExpn;
}
/**
* @param nIndExpn_ The dimensionality of the function returned by {@link #calculateFutureConditionalExpectations(int[], int[], State)}.
*/
final protected void setNumberOfIndividualExpectations(int nIndExpn_)
{
_nIndExpn = nIndExpn_;
}
/**
* Override this method to perform validation checks, e.g. for monotonicity,
* on the calculated implied start-of-period state. Note that if the
* validation is expensive, it would be better to only perform it in some
* periods.
*
* The default implementation does nothing.
*
* @param impliedStartOfPeriodState_
* The calculated implied start-of-period state
* @param state_
* The current state of the calculation
*
* @throws ModelException
* If the validation fails
*/
protected void validateStartOfPeriodState(DoubleArray impliedStartOfPeriodState_, S state_) throws ModelException
{
}
/**
* Override this method to perform validation checks, e.g. for monotonicity,
* on the calculated expected future conditions. Note that if the validation
* is expensive, it would be better to only perform it in some periods.
*
* The default implementation does nothing.
*
* @param expectedFutureCondition_ The calculated expected future condition to validate
* @param state_ The current state of the calculation
*
* @throws ModelException If the validation fails
*/
protected void validateExpectedFutureCondition(DoubleArray expectedFutureCondition_, S state_) throws ModelException
{
}
/**
* Calculates the future individual conditional expectations given the
* current and future exogenous states, both individual and aggregate, and
* the current aggregate variables
*
* @param exogenousStates_ The indexes of current and future aggregate and future individual
* exogenous states to perform the calculation for
* @param aggEndoStateIndex_ The indexes of current aggregate endogenous states to perform the
* calculation for
* @param state_ The calculation state
*
* @return A function yielding the expectational part of the Euler equation(s) given future individual
* endogenous states, conditional on the input variables
*
* @throws ModelException If model-specific processing errors occur
*/
protected abstract DoubleArray calculateFutureConditionalExpectations(int[] exogenousStates_, int[] aggEndoStateIndex_, S state_) throws ModelException;
/**
* Given the expectation of the value of the future part of the inter-temporal condition, conditional on the
* present end-of-period state, calculate the implied present start-of-period state
*
* @param expectedConditions_ The future value of the intertemporal condition, conditional on the present end-of-period state
* @param state_ The state of the current calculation
*
* @return The implied present start-of-period state
* @throws ModelException If model-specific processing errors occur
*/
protected abstract DoubleArray> calculateImpliedStartOfPeriodState(DoubleArray expectedConditions_, S state_) throws ModelException;
}