Add comprehensive solution algorithms documentation page

[Copilot]

This PR adds a new documentation page (docs/src/tutorials/solution_algorithms.md) that provides detailed explanations of the perturbation solution algorithms implemented in MacroModelling.jl for solving DSGE models at first, second, and third order.

Overview

The documentation addresses the need for in-depth technical documentation of the solution methods, including exact mathematical formulations, step-by-step algorithms, and matrix construction details.

Content

Fundamental Equations with Nested Taylor Expansion Structure

Each perturbation order section begins with "The Equation to Solve" that clearly shows the nested structure of the Taylor expansion:

• First-order: The fundamental condition that S₁ must satisfy when substituted into the linearized model equations
• Second-order: Shows the first-order differential operator acting on S₂ plus new second-order terms from the Hessian (∇₂). Includes the equivalent generalized Sylvester equation form A S₂ B + C = S₂
• Third-order: Shows the same first-order operator acting on S₃, plus third-order terms (∇₃), plus second-order interaction terms (∇₂ with S₁ and S₂), plus expectation effects (∇₊ S₂)

Each section includes "Key insights" explaining that lower-order conditions are already satisfied, so higher-order solutions only balance new residual terms. This clearly demonstrates that:

• The same first-order differential operator appears at all orders
• Higher-order solutions build on lower-order solutions (nested structure)
• Each order solves for corrections not captured by lower orders

Tensor Representation and Flattening

Comprehensive explanation of how the algorithm works in flattened 2D space:

• Natural tensor form: How second and third-order solutions are 3-tensors and 4-tensors respectively
• Symmetry and compression: How partial derivative symmetry enables compression (e.g., from n²ₛ to nₛ(nₛ+1)/2 entries for second-order)
• Flattening to 2D matrices: Mapping between tensor and flattened representations using Kronecker products, with worked example for n_s=3
• Compression/uncompression matrices:
  • C₂, U₂: Compression and uncompression for second-order (with explicit 9×6 matrix example)
  • C₃, U₃: Compression and uncompression for third-order with three-way symmetry
• Permutation matrices: Complete documentation of all permutation matrices (P, P₁ₗ, P₁ᵣ, P₁ₗ̂, P₂ₗ̂, P₁ₗ̄, P₂ₗ̄, P₁ᵣ̃, P₂ᵣ̃, SP) and their roles in handling different index orderings

Detailed Algorithm Documentation

First-Order Perturbation Solution

• The equation S₁ must satisfy
• Solving the Quadratic Matrix Equation (QME) for dynamic variables
• QR decomposition for numerical stability
• LU factorization for static variables and shock responses
• Multiple solver algorithms (Schur decomposition, doubling algorithm, iterative methods)
• Matrix construction with explicit dimensions

Second-Order Perturbation Solution

• The equation S₂ must satisfy showing nested structure
• Constructing augmented state vectors with perturbation parameter
• Setting up the generalized Sylvester equation A S₂ B + C = S₂ in compressed space
• Exploiting symmetry through compression matrices
• Sparse storage considerations
• Matrix construction accounting for m₂ = nₛ(nₛ+1)/2 unique entries

Third-Order Perturbation Solution

• The equation S₃ must satisfy with all S₂ interaction terms
• Complex Sylvester equation with permutation matrices in compressed space
• Efficient computation of third-order Kronecker products
• Multiple permutations handling different index orderings
• S₂ interactions: coupled with S₁ through various Kronecker products, multiplied by ∇₊ for expectation effects, appearing in Hessian terms through ∇₂
• Three-way symmetry exploitation
• Matrix construction for m₃ = nₛ(nₛ+1)(nₛ+2)/6 unique entries

Additional Topics

• Pruning: Implementation of the pruned state-space system (Andreasen et al. 2018)
• Computational considerations: Automatic differentiation, sparsity exploitation, numerical stability, caching
• Practical example: Working code demonstrating solution computation at different orders
• Academic references: Citations to foundational papers (Villemot 2011, Andreasen et al. 2018, Levintal 2017, Dynare)

Integration

The page is integrated into docs/make.jl as a tutorial, positioned after basic tutorials and before advanced calibration and estimation tutorials.

Example Usage

using MacroModelling

@model RBC begin
    1/c[0] = (β/c[1]) * (α * exp(z[1]) * k[0]^(α-1) + (1-δ))
    c[0] + k[0] = (1-δ) * k[-1] + q[0]
    q[0] = exp(z[0]) * k[-1]^α
    z[0] = ρᶻ * z[-1] + σᶻ * ϵᶻ[x]
end

@parameters RBC begin
    α = 0.5
    β = 0.95
    δ = 0.02
    ρᶻ = 0.9
    σᶻ = 0.01
end

# Get different order solutions
get_solution(RBC)  # First order (default)
get_solution(RBC, order = 2)  # Second order
get_solution(RBC, order = 3)  # Third order

This documentation provides users with complete understanding from theory to implementation: the fundamental equations with proper nested Taylor expansion structure, how tensor problems are translated to the flattened 2D space (including compression and permutation matrices), and the step-by-step algorithms to compute solutions efficiently.

Original prompt

create a page in the documentation where you go through the solution of the model. write down the exact formulations, algorithms, matrices and their construction. do this for first order, second order and third order

---

💡 You can make Copilot smarter by setting up custom instructions, customizing its development environment and configuring Model Context Protocol (MCP) servers. Learn more Copilot coding agent tips in the docs.

[thorek1]

@copilot the higher order algorithm works in the flattened 2d space. can you lay out how the original problem formulation in terms of tensors looks like and how that is translated to the flattened space (incl what the permutation matrices do)

[Copilot]

@copilot the higher order algorithm works in the flattened 2d space. can you lay out how the original problem formulation in terms of tensors looks like and how that is translated to the flattened spa...

Added a comprehensive section "Tensor Representation and Flattening" that explains:

1. Natural tensor form: How second and third-order solutions are 3-tensors and 4-tensors respectively, with examples showing the indexing

2. Symmetry and compression: Why we can compress tensors (e.g., from n²ₛ to n_s(n_s+1)/2 for second-order) due to partial derivative symmetry

3. Flattening process: How tensors are mapped to 2D matrices using Kronecker products, with concrete example for n_s=3

4. Compression matrices (C₂, C₃): Select unique symmetric entries from full Kronecker products

5. Uncompression matrices (U₂, U₃): Expand compressed solutions back to full form for multiplication

6. Permutation matrices: Detailed explanation of all permutation matrices:

   • P: Combines different index orderings in derivatives
   • P₁ₗ, P₁ᵣ: Basic left/right permutations for index swapping
   • P₁ₗ̂, P₂ₗ̂: Permutations in full derivative space dimension
   • P₁ₗ̄, P₂ₗ̄: Permutations in compressed state space
   • P₁ᵣ̃, P₂ᵣ̃: Right permutations for compressed operations
   • SP: Sparse permutation exploiting derivative sparsity

The section includes a worked example showing the explicit 9×6 compression matrix C₂ for a 3-variable case.

Commit: 5f39819

[thorek1]

@copilot for the first order solution there is an equation the first order solution must fulfil and then there is the algorithm to get to the matrix that fulfils that equation. write that equation down at the beginning of the section explaining the first order algorithm and then do the same for the higher order algorithms

[Copilot]

@copilot for the first order solution there is an equation the first order solution must fulfil and then there is the algorithm to get to the matrix that fulfils that equation. write that equation dow...

Added "The Equation to Solve" subsections at the beginning of each perturbation order section:

First-order: Shows the fundamental condition that S₁ must satisfy when substituted into the linearized model equations

Second-order: Presents both the direct form involving ∇₂ (Hessian) and the equivalent generalized Sylvester equation A S₂ B + C = S₂

Third-order: Shows the equation involving ∇₃ (third-order derivatives), permutations, and coupling with S₂, plus the equivalent Sylvester form

Each section now clearly separates:

1. The mathematical equation the solution must fulfill
2. How it's derived from Taylor expansion
3. The algorithm to actually compute it

This makes it clear that the algorithm sections describe methods to find matrices satisfying these specific equations.

Commit: 3472495

[thorek1]

@copilot shouldn't there be more elements from the second order solution in "The third-order solution matrix S_3 must satisfy"

[Copilot]

@copilot shouldn't there be more elements from the second order solution in "The third-order solution matrix S_3 must satisfy"

Expanded the third-order equation to include additional S₂ interaction terms:

1. Added ∇₁ S₂ (kron(S₁, S₂) + kron(S₂, S₁)) term representing how second-order terms affect future expectations through ∇₊

2. Included "perms" notation in the ∇₂ term to indicate multiple permutations of S₁ and S₂ interactions

3. Enhanced the C matrix description with:

   • ∇₊ S₂ kron(sₜ, s₂) term
   • Multiple ∇₂ interactions: kron(S₁, S₂ kron(S₁, S₁)), kron(S₁, S₂ σ)
   • Definition of s₂ = S₂ kron(s_{t-1}, s_{t-1}) showing second-order state evolution

The equation now reflects that S₂ appears in several forms: coupled with S₁ through Kronecker products (with various permutations), multiplied by the Jacobian ∇₊ to capture expectation effects, and in Hessian terms through ∇₂.

Commit: a57a5fe

[thorek1]

@copilot since this is a Taylor expansion shouldn't all of this be nested, as in I should find the first order condition (equation the solution must satisfy) in the higher order ones...?

[thorek1]

@copilot since this is a Taylor expansion shouldn't all of this be nested, as in I should find the first order condition (equation the solution must satisfy) in the higher order ones...?

[Copilot]

@copilot since this is a Taylor expansion shouldn't all of this be nested, as in I should find the first order condition (equation the solution must satisfy) in the higher order ones...?

Restructured the equations to show the nested Taylor expansion structure:

Second-order: Now explicitly shows the first-order operator (∇₊ S₁ [S₁ʸ, 0]ᵀ + ∇₀) acting on S₂, plus the new second-order terms from ∇₂. Added a "Key insight" note explaining that the first-order condition for S₁ is already satisfied.

Third-order: Shows the same first-order operator acting on S₃, plus third-order terms (∇₃ with S₁), plus second-order interaction terms (∇₂ with S₁ and S₂), plus expectation effects (∇₊ S₂). Added a "Key insight" emphasizing that both first and second-order conditions are already satisfied, so we only balance new third-order terms.

The equations now clearly show:

1. The same first-order differential operator appears in all orders
2. Each higher order adds new terms while lower-order conditions are already satisfied
3. The nested structure of the Taylor expansion

This makes the mathematical structure much clearer - higher-order solutions don't re-solve the lower-order problems, they just add corrections.

Commit: 110c3b1

@copilot since this is a Taylor expansion shouldn't all of this be nested, as in I should find the first order condition (equation the solution must satisfy) in the higher order ones...?

Restructured the equations to show the nested Taylor expansion structure:

Second-order: Now explicitly shows the first-order operator (∇₊ S₁ [S₁ʸ, 0]ᵀ + ∇₀) acting on S₂, plus the new second-order terms from ∇₂. Added a "Key insight" note explaining that the first-order condition for S₁ is already satisfied.

Third-order: Shows the same first-order operator acting on S₃, plus third-order terms (∇₃ with S₁), plus second-order interaction terms (∇₂ with S₁ and S₂), plus expectation effects (∇₊ S₂). Added a "Key insight" emphasizing that both first and second-order conditions are already satisfied, so we only balance new third-order terms.

The equations now clearly show:

1. The same first-order differential operator appears in all orders
2. Each higher order adds new terms while lower-order conditions are already satisfied
3. The nested structure of the Taylor expansion

This makes the mathematical structure much clearer - higher-order solutions don't re-solve the lower-order problems, they just add corrections.

Commit: 110c3b1