diff --git a/CMakeLists.txt b/CMakeLists.txt index 4bedd88..a3e46ed 100644 --- a/CMakeLists.txt +++ b/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} ) diff --git a/include/finmath/StatisticalModels/garch.h b/include/finmath/StatisticalModels/garch.h new file mode 100644 index 0000000..cdd5e57 --- /dev/null +++ b/include/finmath/StatisticalModels/garch.h @@ -0,0 +1,55 @@ +#ifndef GARCH_H +#define GARCH_H + +#include +#include + +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 \ No newline at end of file diff --git a/include/finmath/StatisticalModels/pca.h b/include/finmath/StatisticalModels/pca.h new file mode 100644 index 0000000..697c639 --- /dev/null +++ b/include/finmath/StatisticalModels/pca.h @@ -0,0 +1,48 @@ +#ifndef PCA_H +#define PCA_H + +#include +#include + +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 \ No newline at end of file diff --git a/src/cpp/StatisticalModels/garch.cpp b/src/cpp/StatisticalModels/garch.cpp new file mode 100644 index 0000000..faf862b --- /dev/null +++ b/src/cpp/StatisticalModels/garch.cpp @@ -0,0 +1,127 @@ +#include "finmath/StatisticalModels/garch.h" +#include + +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 \ No newline at end of file diff --git a/src/cpp/StatisticalModels/pca.cpp b/src/cpp/StatisticalModels/pca.cpp new file mode 100644 index 0000000..850c87b --- /dev/null +++ b/src/cpp/StatisticalModels/pca.cpp @@ -0,0 +1,81 @@ +#include "finmath/StatisticalModels/pca.h" +#include + +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_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 \ No newline at end of file diff --git a/test/OptionPricing/binomial_option_pricing_test.cpp b/test/OptionPricing/binomial_option_pricing_test.cpp index 2d98089..aa4c2af 100644 --- a/test/OptionPricing/binomial_option_pricing_test.cpp +++ b/test/OptionPricing/binomial_option_pricing_test.cpp @@ -13,7 +13,115 @@ int binomial_option_pricing_tests() { double tolerance = 0.001; - std::cout << "Binomial-Tree Tests Passed!" << std::endl; + // Test Case 1: Basic European Call Option + { + double S0 = 100.0; // Initial stock price + double K = 100.0; // Strike price + double T = 1.0; // Time to maturity (1 year) + double r = 0.05; // Risk-free rate + double sigma = 0.2; // Volatility + long N = 100; // Number of time steps + + double call_price = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma, N); + std::cout << "Test Case 1 - Basic Call Option Price: " << call_price << std::endl; + assert(call_price > 0.0 && call_price < S0); + } + + // Test Case 2: Basic European Put Option + { + double S0 = 100.0; + double K = 100.0; + double T = 1.0; + double r = 0.05; + double sigma = 0.2; + long N = 100; + + double put_price = binomial_option_pricing(OptionType::PUT, S0, K, T, r, sigma, N); + std::cout << "Test Case 2 - Basic Put Option Price: " << put_price << std::endl; + assert(put_price > 0.0 && put_price < K * std::exp(-r * T)); + } + + // Test Case 3: Deep In-The-Money Call + { + double S0 = 100.0; + double K = 80.0; // Deep ITM + double T = 1.0; + double r = 0.05; + double sigma = 0.2; + long N = 100; + + double call_price = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma, N); + std::cout << "Test Case 3 - Deep ITM Call Price: " << call_price << std::endl; + assert(call_price > (S0 - K * std::exp(-r * T))); + } + + // Test Case 4: Deep Out-of-The-Money Put + { + double S0 = 100.0; + double K = 120.0; // Deep OTM + double T = 1.0; + double r = 0.05; + double sigma = 0.2; + long N = 100; + + double put_price = binomial_option_pricing(OptionType::PUT, S0, K, T, r, sigma, N); + std::cout << "Test Case 4 - Deep OTM Put Price: " << put_price << std::endl; + assert(put_price > 0.0 && put_price < K * std::exp(-r * T)); + } + + // Test Case 5: Convergence with Different Time Steps + { + double S0 = 100.0; + double K = 100.0; + double T = 1.0; + double r = 0.05; + double sigma = 0.2; + + // Test with different number of time steps + long N1 = 50; + long N2 = 100; + long N3 = 200; + + double price1 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma, N1); + double price2 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma, N2); + double price3 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma, N3); + + std::cout << "Test Case 5 - Convergence Test:" << std::endl; + std::cout << "N=50: " << price1 << std::endl; + std::cout << "N=100: " << price2 << std::endl; + std::cout << "N=200: " << price3 << std::endl; + + // Prices should converge as N increases + assert(std::abs(price2 - price1) > std::abs(price3 - price2)); + } + + // Test Case 6: Different Volatility Levels + { + double S0 = 100.0; + double K = 100.0; + double T = 1.0; + double r = 0.05; + long N = 100; + + // Test with different volatility levels + double sigma1 = 0.1; // Low volatility + double sigma2 = 0.2; // Medium volatility + double sigma3 = 0.3; // High volatility + + double price1 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma1, N); + double price2 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma2, N); + double price3 = binomial_option_pricing(OptionType::CALL, S0, K, T, r, sigma3, N); + + std::cout << "Test Case 6 - Volatility Test:" << std::endl; + std::cout << "sigma=0.1: " << price1 << std::endl; + std::cout << "sigma=0.2: " << price2 << std::endl; + std::cout << "sigma=0.3: " << price3 << std::endl; + + // Higher volatility should lead to higher option prices + assert(price1 < price2 && price2 < price3); + } + + std::cout << "All Binomial-Tree Tests Passed!" << std::endl; return 0; } diff --git a/test/StatisticalModels/garch_test.cpp b/test/StatisticalModels/garch_test.cpp new file mode 100644 index 0000000..7d2f012 --- /dev/null +++ b/test/StatisticalModels/garch_test.cpp @@ -0,0 +1,93 @@ +#include +#include +#include +#include "finmath/StatisticalModels/garch.h" + +int garch_tests() { + // Test Case 1: Simple GARCH(1,1) with known parameters + { + // Generate synthetic data with known parameters + int n = 1000; + Eigen::VectorXd returns(n); + double omega = 0.0001; + double alpha = 0.1; + double beta = 0.8; + + // Initialize first return + returns(0) = 0.01; + + // Generate returns using GARCH(1,1) process + for (int t = 1; t < n; ++t) { + double sigma_t = std::sqrt(omega + alpha * returns(t-1) * returns(t-1) + beta * sigma_t * sigma_t); + returns(t) = sigma_t * (rand() / double(RAND_MAX) - 0.5); + } + + // Fit GARCH model + finmath::StatisticalModels::GARCH garch; + garch.fit(returns); + + // Check parameter estimates + double omega_est = garch.get_omega(); + double alpha_est = garch.get_alpha(); + double beta_est = garch.get_beta(); + + // Parameters should be close to true values + assert(std::abs(omega_est - omega) < 0.0001); + assert(std::abs(alpha_est - alpha) < 0.1); + assert(std::abs(beta_est - beta) < 0.1); + + // Check volatility prediction + double vol_pred = garch.predict_volatility(1); + assert(vol_pred > 0); + } + + // Test Case 2: Realistic financial returns + { + // Generate more realistic returns with volatility clustering + int n = 2000; + Eigen::VectorXd returns(n); + + // Initialize with random values + for (int t = 0; t < n; ++t) { + returns(t) = 0.01 * (rand() / double(RAND_MAX) - 0.5); + } + + // Add some volatility clustering + for (int t = 500; t < 1000; ++t) { + returns(t) *= 2.0; // Higher volatility period + } + + // Fit GARCH model + finmath::StatisticalModels::GARCH garch; + garch.fit(returns); + + // Check model diagnostics + Eigen::VectorXd residuals = garch.get_residuals(); + double aic = garch.get_aic(); + double bic = garch.get_bic(); + + // Basic checks + assert(residuals.size() == n); + assert(aic < 0); // AIC should be negative for good fit + assert(bic < aic); // BIC should be more negative than AIC + + // Check volatility predictions + double vol_1 = garch.predict_volatility(1); + double vol_5 = garch.predict_volatility(5); + double vol_10 = garch.predict_volatility(10); + + // Volatility predictions should be positive and decreasing with horizon + assert(vol_1 > 0); + assert(vol_5 > 0); + assert(vol_10 > 0); + assert(vol_1 >= vol_5); + assert(vol_5 >= vol_10); + } + + std::cout << "All GARCH Tests Passed!" << std::endl; + return 0; +} + +int main() { + return garch_tests(); +} \ No newline at end of file diff --git a/test/StatisticalModels/pca_test.cpp b/test/StatisticalModels/pca_test.cpp new file mode 100644 index 0000000..42d3a43 --- /dev/null +++ b/test/StatisticalModels/pca_test.cpp @@ -0,0 +1,54 @@ +#include +#include +#include +#include "finmath/StatisticalModels/pca.h" + +int pca_tests() { + // Test Case 1: Simple 2D data + { + Eigen::MatrixXd data(4, 2); + data << 1, 2, + 2, 4, + 3, 6, + 4, 8; + + finmath::StatisticalModels::PCA pca(data, 1); + + // Check explained variance ratio + Eigen::VectorXd var_ratio = pca.get_explained_variance_ratio(); + assert(var_ratio(0) > 0.99); // Should explain almost all variance + + // Check transformation and inverse transformation + Eigen::MatrixXd transformed = pca.transform(data); + Eigen::MatrixXd reconstructed = pca.inverse_transform(transformed); + + // Check reconstruction error + double error = (data - reconstructed).norm() / data.norm(); + assert(error < 0.01); + } + + // Test Case 2: Random 3D data + { + Eigen::MatrixXd data = Eigen::MatrixXd::Random(100, 3); + finmath::StatisticalModels::PCA pca(data, 2); + + // Check number of components + assert(pca.get_n_components() == 2); + + // Check explained variance ratio sums to less than 1 + Eigen::VectorXd var_ratio = pca.get_explained_variance_ratio(); + assert(var_ratio.sum() < 1.0); + + // Check transformation preserves dimensionality + Eigen::MatrixXd transformed = pca.transform(data); + assert(transformed.cols() == 2); + assert(transformed.rows() == 100); + } + + std::cout << "All PCA Tests Passed!" << std::endl; + return 0; +} + +int main() { + return pca_tests(); +} \ No newline at end of file