pca, garch, and tests for binomial tree

CMakeLists.txt

@@ -12,6 +12,9 @@ set(CMAKE_BUILD_TYPE Release)
 # Add include directories
 include_directories(${PROJECT_SOURCE_DIR}/include)
 
+# Find Eigen package
+find_package(Eigen3 3.3 REQUIRED NO_MODULE)
+
 # Source files
 file(GLOB SOURCES "src/cpp/*/*.cpp")
 
@@ -25,20 +28,31 @@ add_library(finmath_library SHARED ${SOURCES}
     "include/finmath/TimeSeries/simple_moving_average.h"
     "include/finmath/TimeSeries/rsi.h"
     "include/finmath/TimeSeries/ema.h"
+    "include/finmath/StatisticalModels/pca.h"
+    "include/finmath/StatisticalModels/garch.h")
+
+# Link Eigen to the main library
+target_link_libraries(finmath_library PUBLIC Eigen3::Eigen)
     "include/finmath/TimeSeries/rolling_std_dev.h")
 
 # Test executables
 add_executable(black_scholes_test test/OptionPricing/black_scholes_test.cpp)
-target_link_libraries(black_scholes_test finmath_library)
+target_link_libraries(black_scholes_test finmath_library Eigen3::Eigen)
 
 add_executable(binomial_option_pricing_test test/OptionPricing/binomial_option_pricing_test.cpp)
-target_link_libraries(binomial_option_pricing_test finmath_library)
+target_link_libraries(binomial_option_pricing_test finmath_library Eigen3::Eigen)
 
 add_executable(compound_interest_test test/InterestAndAnnuities/compound_interest_test.cpp)
-target_link_libraries(compound_interest_test finmath_library)
+target_link_libraries(compound_interest_test finmath_library Eigen3::Eigen)
 
 add_executable(rsi_test test/TimeSeries/rsi_test.cpp)
-target_link_libraries(rsi_test finmath_library)
+target_link_libraries(rsi_test finmath_library Eigen3::Eigen)
+
+add_executable(pca_test test/StatisticalModels/pca_test.cpp)
+target_link_libraries(pca_test finmath_library Eigen3::Eigen)
+
+add_executable(garch_test test/StatisticalModels/garch_test.cpp)
+target_link_libraries(garch_test finmath_library Eigen3::Eigen)
 
 add_executable(rolling_std_dev_test test/TimeSeries/rolling_std_dev_test.cpp)
 target_link_libraries(rolling_std_dev_test finmath_library)
@@ -54,6 +68,8 @@ add_test(NAME BlackScholesTest COMMAND black_scholes_test)
 add_test(NAME BinomialOptionPricingTest COMMAND binomial_option_pricing_test)
 add_test(NAME CompoundInterestTest COMMAND compound_interest_test)
 add_test(NAME RSITest COMMAND rsi_test)
+add_test(NAME PCATest COMMAND pca_test)
+add_test(NAME GARCHTest COMMAND garch_test)
 add_test(NAME RollingStdDevTest COMMAND rolling_std_dev_test)
 
 # Add a custom target to run all tests
@@ -62,6 +78,8 @@ add_custom_target(build_and_test
     COMMAND ${CMAKE_COMMAND} --build . --target binomial_option_pricing_test
     COMMAND ${CMAKE_COMMAND} --build . --target compound_interest_test
     COMMAND ${CMAKE_COMMAND} --build . --target rsi_test
+    COMMAND ${CMAKE_COMMAND} --build . --target pca_test
+    COMMAND ${CMAKE_COMMAND} --build . --target garch_test
     COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
     WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
 )

include/finmath/StatisticalModels/garch.h

@@ -0,0 +1,55 @@
+#ifndef GARCH_H
+#define GARCH_H
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace finmath {
+namespace StatisticalModels {
+
+class GARCH {
+public:
+    // Constructor with default parameters (GARCH(1,1))
+    GARCH(double omega = 0.0001, double alpha = 0.1, double beta = 0.8);
+
+    // Fit the model to data
+    void fit(const Eigen::VectorXd& returns);
+
+    // Predict volatility
+    double predict_volatility(int steps_ahead = 1) const;
+
+    // Get model parameters
+    double get_omega() const;
+    double get_alpha() const;
+    double get_beta() const;
+
+    // Get fitted values
+    Eigen::VectorXd get_fitted_values() const;
+
+    // Get residuals
+    Eigen::VectorXd get_residuals() const;
+
+    // Get AIC (Akaike Information Criterion)
+    double get_aic() const;
+
+    // Get BIC (Bayesian Information Criterion)
+    double get_bic() const;
+
+private:
+    double omega_;  // Constant term
+    double alpha_;  // ARCH parameter
+    double beta_;   // GARCH parameter
+
+    Eigen::VectorXd fitted_values_;
+    Eigen::VectorXd residuals_;
+    Eigen::VectorXd conditional_variances_;
+
+    // Helper functions
+    void update_parameters(const Eigen::VectorXd& returns);
+    double calculate_likelihood() const;
+};
+
+} // namespace StatisticalModels
+} // namespace finmath
+
+#endif // GARCH_H 

include/finmath/StatisticalModels/pca.h

@@ -0,0 +1,48 @@
+#ifndef PCA_H
+#define PCA_H
+
+#include <vector>
+#include <Eigen/Dense>
+
+namespace finmath {
+namespace StatisticalModels {
+
+class PCA {
+public:
+    PCA(const Eigen::MatrixXd& data, int num_components = 0);
+
+    // Get principal components
+    Eigen::MatrixXd get_components() const;
+
+    // Get explained variance ratio
+    Eigen::VectorXd get_explained_variance_ratio() const;
+
+    // Transform data to principal component space
+    Eigen::MatrixXd transform(const Eigen::MatrixXd& data) const;
+
+    // Inverse transform from principal component space
+    Eigen::MatrixXd inverse_transform(const Eigen::MatrixXd& transformed_data) const;
+
+    // Get number of components
+    int get_n_components() const;
+
+    // Get mean of original data
+    Eigen::VectorXd get_mean() const;
+
+    // Get standard deviation of original data
+    Eigen::VectorXd get_std() const;
+
+private:
+    Eigen::MatrixXd components_;
+    Eigen::VectorXd explained_variance_ratio_;
+    Eigen::VectorXd mean_;
+    Eigen::VectorXd std_;
+    int n_components_;
+
+    void fit(const Eigen::MatrixXd& data);
+};
+
+} // namespace StatisticalModels
+} // namespace finmath
+
+#endif // PCA_H 

src/cpp/StatisticalModels/garch.cpp

@@ -0,0 +1,127 @@
+#include "finmath/StatisticalModels/garch.h"
+#include <cmath>
+
+namespace finmath {
+namespace StatisticalModels {
+
+GARCH::GARCH(double omega, double alpha, double beta)
+    : omega_(omega), alpha_(alpha), beta_(beta) {}
+
+void GARCH::fit(const Eigen::VectorXd& returns) {
+    int n = returns.size();
+
+    // Initialize vectors
+    fitted_values_ = Eigen::VectorXd::Zero(n);
+    residuals_ = Eigen::VectorXd::Zero(n);
+    conditional_variances_ = Eigen::VectorXd::Zero(n);
+
+    // Initialize first conditional variance
+    conditional_variances_(0) = returns.array().square().mean();
+
+    // Update parameters using maximum likelihood
+    update_parameters(returns);
+}
+
+void GARCH::update_parameters(const Eigen::VectorXd& returns) {
+    int n = returns.size();
+
+    // Initialize parameters for optimization
+    Eigen::Vector3d params(omega_, alpha_, beta_);
+
+    // Simple gradient descent optimization
+    double learning_rate = 0.0001;
+    int max_iterations = 1000;
+
+    for (int iter = 0; iter < max_iterations; ++iter) {
+        // Compute conditional variances
+        for (int t = 1; t < n; ++t) {
+            conditional_variances_(t) = omega_ + 
+                alpha_ * returns(t-1) * returns(t-1) + 
+                beta_ * conditional_variances_(t-1);
+        }
+
+        // Compute residuals
+        residuals_ = returns.array() / conditional_variances_.array().sqrt();
+
+        // Compute gradients
+        double grad_omega = 0.0;
+        double grad_alpha = 0.0;
+        double grad_beta = 0.0;
+
+        for (int t = 1; t < n; ++t) {
+            double var_t = conditional_variances_(t);
+            double r_t = returns(t);
+            double r_t_prev = returns(t-1);
+
+            grad_omega += (1.0 / var_t) * (r_t * r_t / var_t - 1.0);
+            grad_alpha += (r_t_prev * r_t_prev / var_t) * (r_t * r_t / var_t - 1.0);
+            grad_beta += (conditional_variances_(t-1) / var_t) * (r_t * r_t / var_t - 1.0);
+        }
+
+        // Update parameters
+        omega_ += learning_rate * grad_omega;
+        alpha_ += learning_rate * grad_alpha;
+        beta_ += learning_rate * grad_beta;
+
+        // Ensure parameters are positive and sum to less than 1
+        omega_ = std::max(0.0, omega_);
+        alpha_ = std::max(0.0, std::min(1.0, alpha_));
+        beta_ = std::max(0.0, std::min(1.0, beta_));
+
+        if (alpha_ + beta_ > 0.99) {
+            double sum = alpha_ + beta_;
+            alpha_ = (alpha_ / sum) * 0.99;
+            beta_ = (beta_ / sum) * 0.99;
+        }
+    }
+}
+
+double GARCH::predict_volatility(int steps_ahead) const {
+    if (steps_ahead <= 0) return std::sqrt(conditional_variances_.tail(1)(0));
+
+    // Compute long-run variance
+    double long_run_var = omega_ / (1.0 - alpha_ - beta_);
+
+    // Compute predicted variance
+    double predicted_var = long_run_var + 
+        std::pow(alpha_ + beta_, steps_ahead) * 
+        (conditional_variances_.tail(1)(0) - long_run_var);
+
+    return std::sqrt(predicted_var);
+}
+
+double GARCH::get_omega() const { return omega_; }
+double GARCH::get_alpha() const { return alpha_; }
+double GARCH::get_beta() const { return beta_; }
+
+Eigen::VectorXd GARCH::get_fitted_values() const {
+    return fitted_values_;
+}
+
+Eigen::VectorXd GARCH::get_residuals() const {
+    return residuals_;
+}
+
+double GARCH::calculate_likelihood() const {
+    int n = residuals_.size();
+    double log_likelihood = -0.5 * n * std::log(2 * M_PI);
+
+    for (int t = 0; t < n; ++t) {
+        log_likelihood -= 0.5 * (std::log(conditional_variances_(t)) + 
+                                residuals_(t) * residuals_(t));
+    }
+
+    return log_likelihood;
+}
+
+double GARCH::get_aic() const {
+    return -2 * calculate_likelihood() + 6; // 3 parameters
+}
+
+double GARCH::get_bic() const {
+    int n = residuals_.size();
+    return -2 * calculate_likelihood() + 3 * std::log(n);
+}
+
+} // namespace StatisticalModels
+} // namespace finmath 

src/cpp/StatisticalModels/pca.cpp

@@ -0,0 +1,81 @@
+#include "finmath/StatisticalModels/pca.h"
+#include <cmath>
+
+namespace finmath {
+namespace StatisticalModels {
+
+PCA::PCA(const Eigen::MatrixXd& data, int num_components) : n_components_(num_components) {
+    fit(data);
+}
+
+void PCA::fit(const Eigen::MatrixXd& data) {
+    // Center the data
+    mean_ = data.colwise().mean();
+    Eigen::MatrixXd centered = data.rowwise() - mean_.transpose();
+
+    // Scale the data
+    std_ = ((centered.adjoint() * centered) / double(data.rows() - 1)).diagonal().array().sqrt();
+    Eigen::MatrixXd scaled = centered.array().rowwise() / std_.transpose().array();
+
+    // Compute covariance matrix
+    Eigen::MatrixXd cov = (scaled.adjoint() * scaled) / double(data.rows() - 1);
+
+    // Compute eigenvalues and eigenvectors
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigen_solver(cov);
+
+    // Sort eigenvalues and eigenvectors in descending order
+    Eigen::VectorXd eigenvalues = eigen_solver.eigenvalues().reverse();
+    Eigen::MatrixXd eigenvectors = eigen_solver.eigenvectors().rowwise().reverse();
+
+    // Set number of components if not specified
+    if (n_components_ == 0) {
+        n_components_ = data.cols();
+    }
+
+    // Store components and explained variance ratio
+    components_ = eigenvectors.leftCols(n_components_);
+    explained_variance_ratio_ = eigenvalues.head(n_components_) / eigenvalues.sum();
+}
+
+Eigen::MatrixXd PCA::get_components() const {
+    return components_;
+}
+
+Eigen::VectorXd PCA::get_explained_variance_ratio() const {
+    return explained_variance_ratio_;
+}
+
+Eigen::MatrixXd PCA::transform(const Eigen::MatrixXd& data) const {
+    // Center and scale the data
+    Eigen::MatrixXd centered = data.rowwise() - mean_.transpose();
+    Eigen::MatrixXd scaled = centered.array().rowwise() / std_.transpose().array();
+
+    // Project onto principal components
+    return scaled * components_;
+}
+
+Eigen::MatrixXd PCA::inverse_transform(const Eigen::MatrixXd& transformed_data) const {
+    // Reconstruct original data
+    Eigen::MatrixXd reconstructed = transformed_data * components_.transpose();
+
+    // Unscale and uncenter
+    reconstructed = reconstructed.array().rowwise() * std_.transpose().array();
+    reconstructed = reconstructed.rowwise() + mean_.transpose();
+
+    return reconstructed;
+}
+
+int PCA::get_n_components() const {
+    return n_components_;
+}
+
+Eigen::VectorXd PCA::get_mean() const {
+    return mean_;
+}
+
+Eigen::VectorXd PCA::get_std() const {
+    return std_;
+}
+
+} // namespace StatisticalModels
+} // namespace finmath