econometron.utils.projection.basis

• ChebyshevBasis Class from econometron.utils.projection.basis

Purpose

The ChebyshevBasis class creates a functional space based on Chebyshev polynomials of the first kind, which are used to approximate solutions to functional equations in DSGE models. It supports operations such as generating Chebyshev nodes, constructing multidimensional grids, evaluating basis functions, computing function approximations, and performing numerical integration via Gauss-Chebyshev and Gauss-Hermite quadrature. The class is designed to work with the ProjectionSolver to solve nonlinear equations by representing policy or value functions as linear combinations of Chebyshev polynomials.

Key Concepts

• Chebyshev Polynomials: Orthogonal polynomials defined on ([-1, 1]) with the recurrence relation: [ T_0(x) = 1, \quad T_1(x) = x, \quad T_n(x) = 2x T_{n-1}(x) - T_{n-2}(x) ] They are ideal for numerical approximation due to their orthogonality and fast convergence properties.
• Role in DSGE Models: Used to approximate policy functions (e.g., consumption, capital) or value functions in DSGE models by representing them as: [ f(x) \approx \sum_{i} c_i \phi_i(x) ] where ( \phi_i(x) ) are multidimensional Chebyshev basis functions, and ( c_i ) are coefficients.
• Projection Methods: The class supports collocation, Galerkin, and least squares methods by providing nodes, basis evaluations, and quadrature for integration.

---

Initialization (__init__(self, order_vector, lower_bounds, upper_bounds))

• Purpose: Initializes the Chebyshev functional space with polynomial orders and domain bounds for each dimension.

• Parameters:

  Name	Type	Description	Default
  order_vector	array-like	Order of Chebyshev polynomials for each dimension.	None
  lower_bounds	array-like	Lower bounds of the domain for each dimension.	None
  upper_bounds	array-like	Upper bounds of the domain for each dimension.	None

• Attributes:

  • self.order_vector (np.ndarray): Polynomial orders (shape: (n_dims,)).
  • self.lower_bounds (np.ndarray): Lower bounds of the domain (shape: (n_dims,)).
  • self.upper_bounds (np.ndarray): Upper bounds of the domain (shape: (n_dims,)).
  • self.n_dims (int): Number of dimensions (length of order_vector).

• Mechanics:

  • Converts input arrays to NumPy arrays with appropriate data types (int for order_vector, float for bounds).
  • Validates that order_vector, lower_bounds, and upper_bounds have the same length.
  • Stores the number of dimensions (n_dims).

• Explanation:

  • Defines the approximation space, where each dimension has a Chebyshev polynomial of order order_vector[i].
  • The domain ([lower_bounds[i], upper_bounds[i]]) specifies the state space (e.g., capital, productivity in a DSGE model).
  • Ensures consistency in dimensionality for subsequent operations.

• Example:

  • For a 2D state space (e.g., capital and productivity) with order_vector=[3, 2], lower_bounds=[0.1, 0.5], upper_bounds=[1.0, 1.5], the class sets up a basis with 4 polynomials (degrees 0 to 3) in the first dimension and 3 (degrees 0 to 2) in the second.

---

Method: funnode(self)

• Purpose: Generates Chebyshev nodes (extrema of Chebyshev polynomials) for collocation methods in each dimension.

• Returns:

  • nodes (list of arrays): List of Chebyshev nodes for each dimension, where each array has shape (order_vector[i] + 1,).

• Mechanics:

  • For each dimension ( i ) with order ( n = \text{order_vector}[i] ):
    • Computes Chebyshev nodes on ([-1, 1]) using: [ x_k = \cos\left(\frac{(2k + 1)\pi}{2(n + 1)}\right), \quad k = 0, 1, \dots, n ]
    • Scales nodes to the domain ([a, b] = [lower_bounds[i], upper_bounds[i]]) using: [ \tilde{x}_k = \frac{b - a}{2} x_k + \frac{a + b}{2} ]
    • Reverses the node order ([::-1]) to match standard ordering (e.g., increasing values).
  • Returns a list of node arrays, one per dimension.

• Mathematical Foundation:

  • Chebyshev nodes are the extrema of the ( n )-th Chebyshev polynomial, optimal for polynomial interpolation due to their distribution, which minimizes Runge’s phenomenon.
  • The transformation maps nodes to the user-specified domain for DSGE state variables.

• Explanation:

  • Used in the collocation method to define points where residuals are set to zero.
  • Ensures high accuracy in polynomial approximations by using optimally spaced nodes.
  • In DSGE models, nodes represent specific states (e.g., capital levels) where equilibrium conditions are enforced.

• Example:

  • For order_vector=[2], lower_bounds=[0], upper_bounds=[1], the nodes are computed as: [ x_k = \cos\left(\frac{(2k + 1)\pi}{6}\right), \quad k = 0, 1, 2 ] Scaled to ([0, 1]), yielding approximately [0.933, 0.5, 0.067] (reversed).

---

Method: gridmake(self, nodes)

• Purpose: Creates a multidimensional grid by taking the Cartesian product of nodes for each dimension.

• Parameters:

  Name	Type	Description
  nodes	list of arrays	List of nodes for each dimension.

• Returns:

  • grid (ndarray): Grid of points with shape (prod(n_i + 1), n_dims), where n_i is the number of nodes in dimension ( i ).

• Mechanics:

  • Uses itertools.product to compute the Cartesian product of the node arrays.
  • Converts the product to a NumPy array, where each row is a point in the multidimensional grid.

• Explanation:

  • The grid is used for evaluating basis functions or residuals at a dense set of points, essential for collocation, Galerkin, or least squares methods.
  • In DSGE models, the grid represents combinations of state variables (e.g., all pairs of capital and productivity levels).
  • The shape ensures all possible combinations of nodes are included.

• Example:

  • For nodes=[[0, 1], [2, 3, 4]], the grid has ( 2 \times 3 = 6 ) points: [ \begin{bmatrix} 0 & 2 \ 0 & 3 \ 0 & 4 \ 1 & 2 \ 1 & 3 \ 1 & 4 \end{bmatrix} ]

---

Method: funbas(self, grid)

• Purpose: Computes the Chebyshev basis functions evaluated at the given grid points.

• Parameters:

  Name	Type	Description
  grid	ndarray	Grid points of shape (n_points, n_dims) or (n_points,) for 1D.

• Returns:

  • basis (ndarray): Basis matrix of shape (n_points, prod(n_i + 1)), where each column corresponds to a multidimensional Chebyshev basis function.

• Mechanics:

  • Ensures grid is 2D; reshapes 1D inputs to (n_points, 1).
  • Validates that grid has n_dims columns.
  • Transforms grid points to ([-1, 1]) for each dimension ( d ): [ x_d = \frac{2 (grid_d - a_d)}{b_d - a_d} - 1, \quad \text{clipped to } [-1, 1] ] where ( a_d, b_d ) are lower_bounds[d], upper_bounds[d].
  • Computes the total basis size: ( \prod (order_vector[i] + 1) ).
  • Initializes basis matrix as ones (shape: (n_points, total_basis_size)).
  • Iterates over all multi-indices (combinations of polynomial degrees) using itertools.product.
  • For each multi-index ( (k_1, k_2, \dots, k_{n_dims}) ), computes the basis function: [ \phi(x) = \prod_{d=1}^{n_dims} T_{k_d}(x_d) ] where ( T_{k_d} ) is the ( k_d )-th Chebyshev polynomial, computed via _chebyshev_poly.
  • Stores the product in the corresponding column of basis.

• Mathematical Foundation:

  • The basis functions are tensor products of Chebyshev polynomials: [ \phi_i(x) = \prod_{d=1}^{n_dims} T_{k_d}(x_d) ] where ( i ) indexes the multi-index ( (k_1, \dots, k_{n_dims}) ).
  • The matrix basis allows linear combinations to approximate functions: ( f(x) \approx \text{basis} \cdot c ).

• Explanation:

  • Essential for evaluating the approximated function or residuals in projection methods.
  • In DSGE models, the basis matrix is used to compute policy function values or weighted residuals at grid points.
  • The transformation to ([-1, 1]) ensures Chebyshev polynomials are evaluated in their natural domain.

• Example:

  • For order_vector=[1, 1], grid=[[0.5, 0.5]], lower_bounds=[0, 0], upper_bounds=[1, 1]:
    • Transform: ( x = [0, 0] ) (since ( (0.5 - 0)/(1 - 0) \cdot 2 - 1 = 0 )).
    • Basis functions: ( T_0(0) = 1, T_1(0) = 0 ), so the basis matrix is: [ \begin{bmatrix} 1 \cdot 1 & 1 \cdot 0 & 0 \cdot 1 & 0 \cdot 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix} ]

---

Method: funeval(self, coefficients, evaluation_points)

• Purpose: Evaluates the approximated function ( f(x) = \sum c_i \phi_i(x) ) at specified points using the Chebyshev coefficients.

• Parameters:

  Name	Type	Description
  coefficients	array-like	Coefficients of the approximation (shape: prod(n_i + 1)).
  evaluation_points	ndarray	Points where the function is evaluated (shape: (n_points, n_dims) or (n_points,) for 1D).

• Returns:

  • values (ndarray): Function values at evaluation points (shape: (n_points,)).

• Mechanics:

  • Reshapes evaluation_points to 2D if 1D.
  • Computes the basis matrix at evaluation_points using self.funbas.
  • Validates that coefficients length matches the basis size.
  • Computes values via matrix multiplication: [ \text{values} = \text{basis} \cdot \text{coefficients} ]

• Explanation:

  • Evaluates the approximated policy or value function at arbitrary points, crucial for simulation or analysis in DSGE models.
  • The basis matrix maps coefficients to function values efficiently.
  • Used by ProjectionSolver.evaluate_solution to compute solutions after solving for coefficients.

• Example:

  • For coefficients=[1, 0, 0, 0], evaluation_points=[[0.5, 0.5]], the output is approximately [1] (based on the basis matrix from the previous example).

---

Method: gauss_chebyshev_quadrature(self, n_nodes)

• Purpose: Computes Gauss-Chebyshev quadrature nodes and weights for numerical integration in each dimension.

• Parameters:

  Name	Type	Description
  n_nodes	int or array-like	Number of quadrature nodes per dimension (single int or list).

• Returns:

  • nodes (list of arrays): Quadrature nodes for each dimension.
  • weights (list of arrays): Quadrature weights for each dimension.

• Mechanics:

  • If n_nodes is an integer, applies it to all dimensions; otherwise, validates length matches n_dims.
  • For each dimension ( i ) with ( n = n_nodes[i] ):
    • Computes nodes on ([-1, 1]): [ x_k = \cos\left(\frac{(2k + 1)\pi}{2n}\right), \quad k = 0, 1, \dots, n-1 ]
    • Computes weights: ( w_k = \frac{\pi}{n} ).
    • Scales nodes to ([a, b]): [ \tilde{x}_k = \frac{b - a}{2} x_k + \frac{a + b}{2} ]
    • Scales weights: ( w_k = w_k \cdot \frac{b - a}{2} ).
    • Reverses node order for consistency.
  • Returns lists of nodes and weights.

• Mathematical Foundation:

  • Gauss-Chebyshev quadrature approximates integrals: [ \int_{-1}^{1} f(x) \frac{1}{\sqrt{1 - x^2}} dx \approx \sum_{k=1}^n w_k f(x_k) ]
  • After scaling to ([a, b]), it approximates: [ \int_a^b f(x) dx ]

• Explanation:

  • Used in the Galerkin method to compute weighted residuals via numerical integration.
  • In DSGE models, quadrature is critical for evaluating expectations or integrals in equilibrium conditions.

• Example:

  • For n_nodes=2, lower_bounds=[0], upper_bounds=[1]:
    • Nodes: ( x_k = \cos(\pi (2k + 1) / 4) ), scaled to [0.8536, 0.1464].
    • Weights: ( w_k = \frac{\pi}{2} \cdot \frac{1 - 0}{2} = \frac{\pi}{4} ).

---

Method: integrate(self, func, n_nodes)

• Purpose: Performs multidimensional Gauss-Chebyshev quadrature to integrate a given function over the domain.

• Parameters:

  Name	Type	Description
  func	callable	Function to integrate, takes array of shape (n_points, n_dims).
  n_nodes	int or array-like	Number of quadrature nodes per dimension.

• Returns:

  • integral (float): Approximate integral value.

• Mechanics:

  • Computes quadrature nodes and weights using self.gauss_chebyshev_quadrature.
  • Creates a grid from nodes using self.gridmake.
  • Evaluates func at grid points.
  • Computes the multidimensional weight as the Kronecker product of per-dimension weights.
  • Calculates the integral: [ \text{integral} = \sum_{\text{grid}} (\text{total_weight} \cdot f(\text{grid})) ]

• Explanation:

  • Approximates integrals like ( \int f(x) dx ) over the multidimensional domain, used in Galerkin methods or DSGE model expectations.
  • The Kronecker product ensures correct weighting for multidimensional integration.

• Example:

  • For func=lambda x: x[:, 0] ** 2, n_nodes=2, the method computes the integral over the domain using quadrature nodes.

---

Method: gauss_hermite_nodes(self, n, sigma)

• Purpose: Computes Gauss-Hermite quadrature nodes and weights for integration over a normal distribution ( N(0, \sigma^2) ), with a transformation to log-normal.

• Parameters:

  Name	Type	Description
  n	int	Number of quadrature nodes.
  sigma	float	Standard deviation of the normal distribution.

• Returns:

  • nodes (ndarray): Quadrature nodes for the log-normal distribution.
  • weights (ndarray): Quadrature weights.

• Mechanics:

  • Uses np.polynomial.hermite.hermgauss to compute nodes and weights for standard Gauss-Hermite quadrature.
  • Scales nodes for ( N(0, \sigma^2) ): [ \text{nodes} = \sqrt{2} \cdot \sigma \cdot \text{nodes} ]
  • Scales weights: [ \text{weights} = \text{weights} / \sqrt{\pi} ]
  • Transforms nodes to log-normal: ( \text{nodes} = \exp(\text{nodes}) ).

• Mathematical Foundation:

  • Gauss-Hermite quadrature approximates integrals: [ \int_{-\infty}^{\infty} f(z) e^{-z^2} dz \approx \sum_{k=1}^n w_k f(z_k) ]
  • Scaled for ( N(0, \sigma^2) ), and transformed for log-normal expectations.

• Explanation:

  • Used in DSGE models to compute expectations over stochastic shocks (e.g., productivity shocks following a log-normal distribution).
  • The log-normal transformation is common for positive-valued shocks.

• Example:

  • For n=3, sigma=1, computes nodes and weights for integrating over ( N(0, 1) ), then transforms to log-normal.

---

Method: _chebyshev_poly(self, x, n)

• Purpose: Computes the ( n )-th Chebyshev polynomial of the first kind at points ( x ).

• Parameters:

  Name	Type	Description
  x	array-like	Points where the polynomial is evaluated.
  n	int	Polynomial degree.

• Returns:

  • poly (ndarray): Polynomial values at x.

• Mechanics:

  • Handles base cases:
    • ( n = 0 ): Returns ones (since ( T_0(x) = 1 )).
    • ( n = 1 ): Returns x (since ( T_1(x) = x )).
  • For ( n \geq 2 ), uses the recurrence relation: [ T_n(x) = 2x T_{n-1}(x) - T_{n-2}(x) ]
  • Iteratively computes ( T_n(x) ) using two previous polynomials.

• Explanation:

  • Core function for evaluating Chebyshev polynomials, used by funbas to construct basis functions.
  • Ensures numerical stability through the recurrence relation.

• Example:

  • For x=[0, 0.5], n=2: [ T_2(x) = 2x^2 - 1 \implies T_2([0, 0.5]) = [-1, -0.5] ]

---

Notes

• Integration with ProjectionSolver:

  • The ChebyshevBasis class provides the mathematical foundation for ProjectionSolver’s collocation, Galerkin, and least squares methods.
  • Methods like funnode, gridmake, and funbas support collocation and Galerkin by providing nodes and basis evaluations.
  • gauss_chebyshev_quadrature and integrate are critical for Galerkin’s weighted residual calculations.
  • gauss_hermite_nodes supports stochastic DSGE models with normal or log-normal shocks.

• Numerical Stability:

  • Chebyshev polynomials are orthogonal, reducing numerical errors in high-degree approximations.
  • Clipping in funbas ensures robustness against points slightly outside the domain.

• Limitations:

  • Curse of Dimensionality: The grid size grows exponentially with n_dims, limiting scalability.
  • Domain Restriction: Assumes a rectangular domain defined by lower_bounds and upper_bounds.
  • Log-Normal Specifics: gauss_hermite_nodes is tailored for log-normal shocks, which may not suit all DSGE models.

• Dependencies:

  • numpy: For array operations and polynomial computations.
  • itertools.product: For generating multidimensional grids and basis indices.
  • np.polynomial.hermite.hermgauss: For Gauss-Hermite quadrature.

---

Example Use Case

import numpy as np
from econometron.utils.projection.basis import ChebyshevBasis

# Initialize Chebyshev basis
order_vector = [3, 2]
lower_bounds = [0.1, 0.5]
upper_bounds = [1.0, 1.5]
basis = ChebyshevBasis(order_vector, lower_bounds, upper_bounds)

# Generate Chebyshev nodes
nodes = basis.funnode()
print("Nodes:", nodes)

# Create grid
grid = basis.gridmake(nodes)
print("Grid shape:", grid.shape)

# Evaluate basis functions
basis_matrix = basis.funbas(grid)
print("Basis matrix shape:", basis_matrix.shape)

# Evaluate function with coefficients
coeffs = np.ones(np.prod(np.array(order_vector) + 1))
values = basis.funeval(coeffs, grid)
print("Function values:", values)

# Perform quadrature
def func(x): return x[:, 0] ** 2 * x[:, 1]
integral = basis.integrate(func, n_nodes=[5, 5])
print("Integral:", integral)

# Gauss-Hermite quadrature
nodes, weights = basis.gauss_hermite_nodes(n=5, sigma=1.0)
print("Log-normal nodes:", nodes)