API Reference

PE_Unit(b_::Float64, base_::Float64, d_::Int)
PE_Unit(b_::Float64, base_::Date, d_::Int)

This creates a PEUnit which has a functional form of exp(b(x-base)) (x-base)^d. They cannot be used in any productive way by themselves but are needed to construct a PEFunction. An empty PEUnit (which might be used to create a constant PEFunction) can be created by PEUnit().

PE_Function(multiplier_::Float64, units_::UnitMap)

This is the main constructor for a PE Function. The functional form of the function is the multiplier multiplied by all PE_Units.

For instance the PEFunction created by PEFunction(6.0, makeunitmap([:x, :y] .=> [PEUnit(1.0,1.0,1), PEUnit(0.0,2.0,4)])) has a functional form of 6 (x-1) exp(x-1) (y-2)^4

The following convenience functions create a PEFunction where there is only one variable (with a symbol :default). PEFunction(multiplier::Float64,b::Float64, base::Float64, d::Int) PEFunction(multiplier::Float64,b::Float64, base::Date, d::Int) The following convenience function creates a PEFunction where there are no variables and hence it is constant: PE_Function(num::Float64 = 0.0)

Sum_Of_Functions(functions)

Creates a SumOfFunctions from an array of PEFunctions and/or SumOfFunctions. The constructors for this type go through each input SumOfFunctions and takes out the contained PEFunctions (so unecessary nesting doesnt occur where a SumOfFunctions could contain another SumOfFunctions). The constructors also aggregate PEFunctions where possible. For intance if two PEFunctions have the same PEUnits and differ in their multiplier these multipliers can be added. The constructors also remove zero multiplier PEFunctions.

Piecewise_Function(functions_::Array{Union{Missing,Sum_Of_Functions}}, thresholds_::OrderedDict{Symbol,Array{Float64,1}})

Creates a PiecewiseFunction from a multidimensional array of SumOfFunctions and an ordered dict of thresholds. The xth dimension in the thresholds dict corresponds to the xth dimension of the array of functions. A function can be lookuped up considering the thresholds and then selecting from the array. For instance if the first dimension is denoted :y and it's thresholds are [-4.0,0.0,3.4] and we query at a point with a :y coordinate of 2.7 then the function we look up will be from the file functions[2, ...] where ... represents the coordinates the the other dimensions. This is because 2.7 is greater than the second element of [-4.0,0.0,3.4] but less than the third. If this piecewise function were to be queried at a :y coordinate of -5.0 then a missing value will be returned. To specify piecewise functions on an unlimited domain the first element of the threshold can be set as -Inf. To set a limited domain on the upper end then add a Missing value to the functions_ array. Any other (ie interior) point can also be made undefined by putting a missing value into the functions_ array.

Note that PiecewiseFunction works by assigning a SumOfFunctions to every region within the space defined by the thresholds dict. It is only possible to specify a region as a hypercube however. More complex regions are not possible.

Note too that piecewise functions will scale poorly in high dimensions. If there are 10 dimensions and each has 4 elements in its threshold dict then the array for the piecewise function will have more than one million entries. In cases where there are no interactions between dimensions it is more efficient to use a SumOfPiecewiseFunctions object (which is basically an array of Piecewise Functions). For instance consider the following function: f(x,y,z) = max(x,5) + max(y,3) + max(z,3) We could code this as a piecewise function or as the sum of three piecewise functions. The three piecewise function implementation will contain fewer PEFunctions.

Sum_Of_Piecewise_Functions(functions_::Array{Piecewise_Function,1}, global_funcs_::Sum_Of_Functions)

At the cost of being less flexible a SumOfPiecewise_Functions is more efficient than a PiecewiseFunction. Use this if trying to represent a piecewise function that can be decomposed into a sum of lower dimensional piecewise functions.

years_between(a::Date, b::Date)
years_between(a::Dates.Day, b::Dates.Day)

Returns the number of years that have elapsed since 1-Jan-2000. For the purposes of this calculation there are 365.2422 days in a year.

years_between(until::Date, from::Date)
years_between(until::Dates.Day, from::Dates.Day)

Returns the number of years between two dates. For the purposes of this calculation there are 365.2422 days in a year.

Period length is designed to convert TimePeriod objects to a float in a consistent way to years_from_global_base

evaluate(f::MultivariateFunction, coordinates::Dict{Symbol,Float64})
evaluates a function at coordinates.

For univariate functions with a variable name of :default (such as those created by PE_Function's convenience functions)
evaluation can take place with no dictionary:
evaluate(f::MultivariateFunction, coordinate::Float64)

evaluate(hess::Hessian, coordinates::Dict{Symbol,Float64})

This evaluates a Hessian object to create a Hermitian matrix representing the hessian at a point.

underlying_dimensions(f::MultivariateFunction)
Returns a set containing all of the variables upon which f depends.

≂

This tests whether two structs are close after allowing for numerical tolerance.

rebadge(f::PE_Function, mapping::Dict{Symbol,Symbol})
can be used to change the names of the variables in a MultivariateFunction.

change_base(f::PE_Function, new_bases::BaseMap)

This function changes the bases in the PEUnits of a PEFunction. This is useful for getting two PEFunctions comformable for simpler multiplication. Often a base change means that an array of PEFunctions are needed to represent a function. So an Array{PE_Function,1} is returned.

+(f::MultivariateFunction,number::Float64) +(f::MultivariateFunction,number::Int)

A Multivariate Function can be added to an scalar to from a new MultivariateFunction. This action promotes a PEFunction to a SumOf_Functions. The type of all other MultivariateFunctions is unchanged.

+(f1::MultivariateFunction, f2::MultivariateFunction)
Any two multivariate Functions can be added to form a MultivariateFunction reflecting the sum.

-(f::MultivariateFunction,number::Float64) -(f::MultivariateFunction,number::Int)

A scalar can be subtracted from a Multivariate Function. This action promotes a PEFunction to a SumOf_Functions. The type of all other MultivariateFunctions is unchanged.

-(f1::MultivariateFunction, f2::MultivariateFunction)
Any MultivariateFunction can be subtracted from another to form a MultivariateFunction reflecting the difference.

*(f::MultivariateFunction,number::Float64) *(f::MultivariateFunction,number::Int)

A Multivariate can be multiplied by a scalar. This does not change the type of any MultivariateFunction.

*(f1::MultivariateFunction, f2::MultivariateFunction)
Any two MultivariateFunctions can be multiplied to form a MultivariateFunction reflecting the product.

/(f::MultivariateFunction,number::Float64) /(f::MultivariateFunction,number::Int)

A Multivariate can be divided by a scalar. This does not change the type of any MultivariateFunction. Note that the opposite operation cannot be done. While f / 5 is permitted 5 / f is not supported by this package. It is not possible to divide by a function.

^(f::MultivariateFunction,number::Int)

A Multivariate can be raised by a positive Integer power. This will generally promote a PEFunction to a SumOf_Functions. Note that the opposite operation cannot be done. While f ^ 5 is permitted 5 ^ f is not supported by this package. It is not possible to raise by the power of a function.

 derivative(f::MultivariateFunction, derivs::Dict{Symbol,Int})

This generates a function representing the derivative of function f. The derivative is that specified by the derivs dict. So if derivs is Dict{[:x,:y] .=> [1,2]} then there will be one derivative with respect to x and 2 with respect to y.

 all_derivatives(f::MultivariateFunction, degree::Int = 2, dimensions::Set{Symbol} = underlying_dimensions(f))

This generates a dict containing functions representing all of the deriviates of a function up to the order of degree.

 integral(f::PE_Function, limits::Dict{Symbol,Tuple{Union{Symbol,Float64},Union{Symbol,Float64}}})
 integral(f::Sum_Of_Functions, limits::Dict{Symbol,Tuple{Union{Symbol,Float64},Union{Symbol,Float64}}})

This gives a function representing the integral of a function, f, with limits in each dimension given by a dict. The dict should contain a tuple for each variable in the function. The left member of the tuple contains the lower limit and the right member the upper limite. Each can be a Float64 or a symbol. If a symbol is input then this will get incorporated as a dimension in the MultivariateFunction created by the integral function.

Hessian(f::MultivariateFunction, dimensions::Array{Symbol,1})
Hessian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels_::Array{Symbol,1})

This represents expressions for constructing a hessian matrix for a function. It can be evaluated to get a Hermitian matrix of the hessian at a particular location.

jacobian(derivs::Union{Dict{Dict{Symbol,Int},PE_Function},Dict{Dict{Symbol,Int},Sum_Of_Functions},Dict{Dict{Symbol,Int},Piecewise_Function}}, labels::Array{Symbol,1})
jacobian(f::MultivariateFunction, dimensions::Array{Symbol,1})

This generates an array of MultivariateFunctions representing the derivatives of a function. This array can be evaluated with evaluate.(jacobian, Ref(coordinates)) to give a vector of the derivative values at a point.

uniroot(f::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)

This takes the analytical jacobian and hessian of a function and uses them to find a nearby root. It finds a root using Newton's method.

find_local_optima(func::MultivariateFunction, initial_guess::Dict{Symbol,Float64}; step_size::Float64 = 1.0, max_iters::Int = 40, convergence_tol::Float64 = 1e-10, print_reports::Bool = false)

This takes the analytical jacobian and hessian of a function and uses them to find a nearby optima. The optima it will find are based on Newton's method. There is no way to specify whether a minimum or a maximum is sought in Newton's method (at least the pure version of it) and thus this function cannot selectively search for a maximum or minimum. It simply searches for a stationary point.

 integral(f::Piecewise_Function, limits::Dict{Symbol,Tuple{Float64,Float64}})
 integral(f::Sum_Of_Piecewise_Functions, limits::Dict{Symbol,Tuple{Float64,Float64}})

This gives a function representing the integral of a function, f, with limits in each dimension given by a dict. The dict should contain a tuple for each variable in the function. The left member of the tuple contains the lower limit and the right member the upper limite. Each must be a Float64 (Support for inputting a symbol is planned but not yet implemented).

create_quadratic_spline(x::Array{Date,1},y::Array{Float64,1} ; gradients::Union{missing,Array{Float64,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve, Curve), left_gradient::Union{Missing,Float64} = missing, right_gradient::Union{Missing,Float64} = missing, dim_name::Symbol = default_symbol)
create_quadratic_spline(x::Array{Int,1},y::Array{Float64,1} ; gradients::Union{missing,Array{Float64,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve, Curve), left_gradient::Union{Missing,Float64} = missing, right_gradient::Union{Missing,Float64} = missing, dim_name::Symbol = default_symbol)
create_quadratic_spline(x::Array{Float64,1},y::Array{Float64,1} ; gradients::Union{missing,Array{Float64,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve, Curve), left_gradient::Union{Missing,Float64} = missing, right_gradient::Union{Missing,Float64} = missing, dim_name::Symbol = default_symbol)
create_quadratic_spline(schum::Schumaker; dim_name::Symbol = default_symbol)

Create a quadratic spline. The spline is a Schumaker shape-preserving spline which is taken from the SchumakerSpline.jl package.

create_linear_interpolation(x::Array{Date,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_linear_interpolation(x::Array{<:DatePeriod,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_linear_interpolation(x::Array{Float64,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)

Create a piecewise linear one-dimensional function which interpolates linearly between datapoints.

create_constant_interpolation_to_left(x::Array{Date,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_constant_interpolation_to_left(x::Array{<:DatePeriod,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_constant_interpolation_to_left(x::Array{Float64,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)

Create a piecewise constant one-dimensional function which carries values from the right to the left.

create_constant_interpolation_to_right(x::Array{Date,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_constant_interpolation_to_right(x::Array{<:DatePeriod,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)
create_constant_interpolation_to_right(x::Array{Float64,1},y::Array{Float64,1}; dim_name::Symbol = default_symbol)

Create a piecewise constant one-dimensional function which carries values from the left to the right.

get_chebyshevs_up_to(num::Int, first_kind::Bool = true; dim_name::Symbol = default_symbol)
get_chevyshevs_up_to(num::Int, first_kind::Bool = true; dim_name::Symbol = default_symbol)
Output all chebyshev polynomials up to degree num.

create_chebyshev_approximation(f::Function, nodes::Int, degree::Int, limits::OrderedDict{Symbol,Tuple{Float64,Float64}}, function_takes_Dict::Bool = false)

Creates a SumOfFunctions that approximates a function, f, with a set of chebyshevs of a particular degree. The nodes input specifies at how many locations the function is to be evaluated for approximation purposes in each dimension. The limits OrderedDict specifies the domain of where the function is to be approximated.

If functiontakesDict is true then the function will be evaluated by inputting a Dict{Symbol,Float64}. Otherwise the function will be evaluated with f(values(limits)...) Note that the order of the OrderedDict specifies the order of inputs to the function in this case.

create_ols_approximation(dd::DataFrame, y::Symbol, model::Array; dropcollinear = true, weights = nothing)
create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Functions; weights = nothing)
create_ols_approximation(dd::DataFrame, y::Symbol, model::Sum_Of_Piecewise_Functions; weights = nothing)

Creates a MultivariateFunction from an OLS regression predicting some variable. You input a dataframe and specify what column in that dataframe is to be predicted by inputting a symbol y. You also input the regression model as an array of basis functions. Each function that is input will be multiplied by the OLS coefficient and will return a new function with these coefficients incorporated.

Returns a tuple (model, regression) where model is a Sum_Of_Piecewise_Functions and regression is the GLM LinearModel object.

If weights is provided (a Vector{Float64} with one non-negative entry per row), a weighted least squares regression is performed using frequency weights.

create_ols_approximation(y::Array{Float64,1}, x::Array{Float64,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_x::Float64 = 0.0, weights = nothing)
create_ols_approximation(y::Array{Float64,1}, x::Array{Date,1}, degree::Int; intercept::Bool = true, dim_name::Symbol = default_symbol, base_date::Date = global_base_date, weights = nothing)

This predicts a linear relationship between the y and x arrays and creates a MultivariateFunction containing the approximation function. The degree specifies how many higher order terms of x should be used (for instance degree 2 implies x and x^2 are both used to predict y).

If weights is provided, a weighted least squares regression is performed.

create_saturated_ols_approximation(dd::DataFrame, y::Symbol, x_variables::Array{Symbol,1}, degree::Int; intercept::Bool = true, bases::Dict{Symbol,Float64} = ..., weights::Union{Nothing, Vector{Float64}} = nothing)

This creates a MultivariateFunction from an OLS regression predicting some variable. You input a dataframe and specify what column in that dataframe is to be predicted by inputting a symbol y. You also put in an array of what x_variables should be used in prediction. A saturated OLS model is then calculated up to the specified degree.

Returns a tuple (model, regression) where model is a Sum_Of_Piecewise_Functions and regression is the GLM LinearModel object.

If weights is provided (a Vector{Float64} with one non-negative entry per row), a weighted least squares regression is performed using frequency weights.

create_mars_spline(dd::DataFrame, y::Symbol, x_variables::Set{Symbol}, MaxM::Int; rel_tol::Float64 = 1e-2, weights = nothing)

This creates a mars spline given a dataframe, response variable and a set of xvariables from the dataframe. The relative tolerance is used in a one-dimensional optimisation step to determine what points at which split values to place a max(0,x-split) function in a particular dimension. The default is intentionally set high because precision is generally not the not that important. For small scale data however you might want to decrease it and increase it for large scale data. You might also want to decrease it if spline creation time doesnt matter much. Note that a small reltol only affects creation time for the spline and not the evaluation time.

Returns a named tuple (model, regression) where model is a Sum_Of_Piecewise_Functions and regression is the GLM LinearModel object.

If weights is provided (a Vector{Float64} with one non-negative entry per row), a weighted least squares fit is used at each step.

create_monotonic_mars_spline(dd::DataFrame, y::Symbol, x_variables::Set{Symbol}, MaxM::Int;
                             rel_tol::Float64 = 1e-2,
                             directions::Dict{Symbol,Symbol} = Dict{Symbol,Symbol}(),
                             min_gradient::Float64 = 0.0,
                             weights::Union{Nothing, Vector{Float64}} = nothing)

Creates a MARS spline that is guaranteed to be monotonic in each specified dimension. Each basis function uses only forward hinges max(0, x - t) (for increasing dimensions) or backward hinges max(0, t - x) (for decreasing dimensions), and non-intercept coefficients are constrained to be non-negative via NNLS fitting. This guarantees the resulting function is monotonic in each dimension in the specified direction.

The directions argument maps each dimension symbol to either :increasing or :decreasing. If not specified, all dimensions default to :increasing.

MaxM specifies the total number of basis functions (including the constant intercept). Each forward step adds one basis function, so MaxM - 1 splits are performed.

If min_gradient is set to a positive value, a linear term with that slope is added in every dimension (with appropriate sign for the direction). This ensures the function is strictly increasing (or decreasing) everywhere with at least the specified gradient, eliminating flat regions. The MARS basis functions are fit to the residual after removing these linear terms.

Returns a named tuple (model, coefficients, rss) where model is a Sum_Of_Piecewise_Functions, coefficients is the NNLS coefficient vector, and rss is the residual sum of squares.

If weights is provided (a Vector{Float64} with one non-negative entry per row), a weighted NNLS fit is used at each step.

create_recursive_partitioning(dd::DataFrame, y::Symbol, x_variables::Set{Symbol}, MaxM::Int; rel_tol::Float64 = 1e-2, weights = nothing)

This creates a recusive partitioning approximation. This seperates the space in to a series of hypercubes each of which has a constant value within the hypercube. Each step of the algorithm divides a hypercube along some dimension so that the different parts of the hypercube can recieve a different value. The relative tolerance is used in a one-dimensional optimisation step to determine what points at which split values to place a hypercube in a particular dimension. The default is intentionally set high because it generally doesnt matter that much. For small scale data however you might want to decrease it and increase it for large scale data. You might also want to decrease it if spline creation time doesnt matter much. Note that a small rel_tol only affects creation time for the spline and not the evaluation time.

Returns a named tuple (model, regression) where model is a Sum_Of_Piecewise_Functions and regression is the GLM LinearModel object.

If weights is provided (a Vector{Float64} with one non-negative entry per row), a weighted least squares fit is used at each step.

trim_mars_spline(dd::DataFrame, y::Symbol, model::Sum_Of_Piecewise_Functions;
                   maximum_RSS::Float64 = -1.0, maximum_increase_in_RSS::Float64 = -1.0,
                   final_number_of_functions::Int = -1, weights = nothing)

This trims a mars spline created in the createmarsspline function. This algorithm goes through each piecewise function in the mars spline and deletes the one that contributes least to the fit. A termination criterion must be set. There are three possible termination criterions. The first is the maximumRSS that can be tolerated. If this is set then functions will be deleted until the deletion of an additional function would push RSS over this amount. The second is maximumincreaseinRSS which will delete functions until a deletion increases RSS by more than this amount. The final is finalnumberof_functions which reduces the number of fucntions to this number.

Returns a named tuple (model, regression) where model is the trimmed Sum_Of_Piecewise_Functions and regression is the final GLM LinearModel object.

If weights is provided (a Vector{Float64} with one non-negative entry per row), weighted RSS is used for all comparisons.

MultivariateFitter

A mutable struct for iteratively fitting multivariate functions. Each call to fit! fits new data and blends with the previous fit. Periodic simplification via trim_mars_spline (on synthetic data from the blended model) reduces complexity.

Fields

• fun - The current fitted MultivariateFunction (or nothing before first fit).
• method - Fitting method. One of :mars, :recursive_partitioning, :monotonic_mars, :ols, :saturated_ols.
• x_variables - Set{Symbol} of predictor dimension names.
• times_through - Number of times fit! has been called.
• MaxM - Number of basis functions for each daily fit.
• simplify_to - Target number of basis functions after trimming. Can be larger than MaxM to allow accumulated complexity. Set to 0 to disable simplification.
• simplification_frequency - Simplify every this many calls to fit!. 0 disables.
• weight_on_new - Maximum blending weight for new fit in (0,1].
• rel_tol - Relative tolerance for MARS/RP optimisation.
• ols_degree - Degree for :saturated_ols method (ignored otherwise).
• directions - Dict{Symbol,Symbol} for monotonic MARS directions.
• min_gradient - Minimum gradient for monotonic MARS (strict monotonicity).

MultivariateAdjustedFitter

A mutable struct for iteratively fitting a shared multivariate shape function with per-group affine adjustments. For each group g, predictions are a_g + b_g * f(x).

Fields

• fun - The current shared MultivariateFunction (or nothing before first fit).
• coefficients - Dict mapping group keys to (a, b) tuples.
• method - Fitting method (same options as MultivariateFitter).
• x_variables - Set{Symbol} of predictor dimension names.
• times_through - Number of times fit! has been called.
• MaxM - Number of basis functions per daily fit.
• simplify_to - Target basis functions after trimming (0 = disabled).
• simplification_frequency - Simplify every N fits (0 = disabled).
• weight_on_new - Maximum blending weight for new fit.
• rel_tol - Optimisation tolerance.
• ols_degree - Degree for :saturated_ols.
• directions - Per-dimension monotonicity directions.
• min_gradient - Minimum gradient for monotonic MARS.
• adjust_for_groups - If true, undo group coefficients before fitting the shape.
• fit_intercept - If true, estimate both a and b; if false, force a=0.
• coefficient_bounds - ((a_min, a_max), (b_min, b_max)) clamping bounds.

fit!(fitter::MultivariateFitter, dd::DataFrame, y::Symbol; weights = nothing)

Fit the MultivariateFitter to new data. The new fit is blended with the previous fit using weight = min(1/(times_through+1), weight_on_new). If simplification_frequency > 0 and the iteration count is a multiple, the model is simplified via backward deletion on synthetic data.

Arguments

• dd - DataFrame containing predictor columns and the response.
• y - Symbol naming the response column.
• weights - Optional Vector{Float64} of non-negative observation weights (one per row). Interpreted as frequency weights. nothing (default) gives uniform weights.

Returns nothing. Access the fitted model via fitter.fun.

fit!(fitter::MultivariateAdjustedFitter, dd::DataFrame, y::Symbol, groups::Vector; weights = nothing)

Fit the MultivariateAdjustedFitter to new data with group labels.

1. Onboard new groups with default coefficients (0.0, 1.0).
2. If adjust_for_groups, undo group coefficients: y_adj = (y - a_g) / b_g.
3. Fit the shared shape function to adjusted data.
4. Blend with previous fit.
5. Re-estimate per-group (a_g, b_g) via weighted OLS, clamped to bounds.
6. Periodically simplify.

Arguments

• dd - DataFrame containing predictor columns and the response.
• y - Symbol naming the response column.
• groups - Vector of group labels (one per row), can be any type (e.g. Symbol, String).
• weights - Optional Vector{Float64} of non-negative observation weights (one per row). Interpreted as frequency weights. nothing (default) gives uniform weights. Used for shape fitting, simplification, and per-group coefficient estimation.

Returns nothing. Access the fitted model via fitter.fun and group coefficients via fitter.coefficients. Groups with fewer than 2 observations are skipped during coefficient estimation.