From a1aeafc84a375243befd6839621fa6a392e33ec1 Mon Sep 17 00:00:00 2001
From: Akshat Sharma <111536616+Akshat111111@users.noreply.github.com>
Date: Tue, 12 Mar 2024 14:35:04 +0530
Subject: [PATCH] Add Monte Carlo simulation, Markowitz portfolio optimization,
 and Kalman filter (#123)

---
 quantitative_finance/kalman_filter.r          | 32 +++++++++++++++++++
 .../markowitz_portfolio_optimization.r        | 28 ++++++++++++++++
 quantitative_finance/monte_carlo_simulation.r | 31 ++++++++++++++++++
 3 files changed, 91 insertions(+)
 create mode 100644 quantitative_finance/kalman_filter.r
 create mode 100644 quantitative_finance/markowitz_portfolio_optimization.r
 create mode 100644 quantitative_finance/monte_carlo_simulation.r

diff --git a/quantitative_finance/kalman_filter.r b/quantitative_finance/kalman_filter.r
new file mode 100644
index 0000000..2967c13
--- /dev/null
+++ b/quantitative_finance/kalman_filter.r
@@ -0,0 +1,32 @@
+library(Metrics)
+set.seed(123)
+num_obs <- 100
+true_returns <- rnorm(num_obs, mean = 0.001, sd = 0.01)
+observed_prices <- cumprod(1 + true_returns) * 100
+noise <- rnorm(num_obs, mean = 0, sd = 0.1)
+noisy_prices <- observed_prices + noise
+# Kalman filter implementation
+kalman_filter <- function(observed_prices) {
+  state <- c(observed_prices[1], 0)
+  P <- matrix(c(1, 0, 0, 1), nrow = 2)
+  Q <- matrix(c(0.0001, 0, 0, 0.0001), nrow = 2)
+  R <- 0.1
+  A <- matrix(c(1, 1, 0, 1), nrow = 2)
+  H <- matrix(c(1, 0), nrow = 1)
+  filtered_states <- matrix(0, nrow = length(observed_prices), ncol = 2)
+  for (i in 1:length(observed_prices)) {
+    state_pred <- A %*% state
+    P_pred <- A %*% P %*% t(A) + Q
+    K <- P_pred %*% t(H) %*% solve(H %*% P_pred %*% t(H) + R)
+    state <- state_pred + K %*% (observed_prices[i] - H %*% state_pred)
+    P <- (matrix(1, nrow = 2, ncol = 2) - K %*% H) %*% P_pred
+    filtered_states[i, ] <- state
+  }
+  return(list(filtered_states = filtered_states, state_pred = state_pred, P_pred = P_pred))
+}
+result <- kalman_filter(noisy_prices)
+plot(observed_prices, type = "l", col = "blue", lwd = 2, main = "Kalman Filter")
+lines(result$filtered_states[, 1], type = "l", col = "red", lwd = 2)
+lines(true_returns, type = "l", col = "green", lwd = 2)
+legend("topright", legend = c("Observed Prices", "Filtered Prices", "True Returns"),
+       col = c("blue", "red", "green"), lty = 1, lwd = 2)
diff --git a/quantitative_finance/markowitz_portfolio_optimization.r b/quantitative_finance/markowitz_portfolio_optimization.r
new file mode 100644
index 0000000..8d831ff
--- /dev/null
+++ b/quantitative_finance/markowitz_portfolio_optimization.r
@@ -0,0 +1,28 @@
+# Required libraries
+library(tidyquant)
+library(quadprog)
+# Set a seed for reproducibility
+set.seed(123)
+# Generate random data for three assets
+num_assets <- 3
+num_obs <- 100
+returns <- matrix(rnorm(num_assets * num_obs), ncol = num_assets)
+# Define the objective function for portfolio optimization
+objective_function <- function(weights, cov_matrix) {
+  portfolio_return <- sum(weights * colMeans(returns))
+  portfolio_volatility <- sqrt(t(weights) %*% cov_matrix %*% weights)
+  return(c(portfolio_return, portfolio_volatility))
+}
+cov_matrix <- cov(returns)
+constraints <- matrix(0, nrow = 2, ncol = num_assets)
+constraints[1, ] <- colMeans(returns)
+constraints[2, ] <- 1
+optimal_weights <- solve.QP(Dmat = 2 * cov_matrix,
+                            dvec = rep(0, num_assets),
+                            Amat = t(constraints),
+                            bvec = c(0.05, 1),
+                            meq = 1)$solution
+cat("Optimal Weights:", optimal_weights, "\n")
+optimal_portfolio <- objective_function(optimal_weights, cov_matrix)
+cat("Expected Return:", optimal_portfolio[1], "\n")
+cat("Volatility:", optimal_portfolio[2], "\n")
diff --git a/quantitative_finance/monte_carlo_simulation.r b/quantitative_finance/monte_carlo_simulation.r
new file mode 100644
index 0000000..de2a819
--- /dev/null
+++ b/quantitative_finance/monte_carlo_simulation.r
@@ -0,0 +1,31 @@
+# Required libraries
+library("quantmod")
+# Parameters
+S0 <- 100    # Initial stock price
+K <- 100     # Strike price
+r <- 0.05    # Risk-free rate
+sigma <- 0.2 # Volatility
+T <- 1       # Time to maturity (in years)
+n <- 252     # Number of trading days
+# Function to simulate stock prices using geometric Brownian motion
+simulate_stock_prices <- function(S0, r, sigma, T, n) {
+  dt <- T/n
+  t <- seq(0, T, by = dt)
+  W <- c(0, cumsum(sqrt(dt) * rnorm(n)))
+  S <- S0 * exp((r - 0.5 * sigma^2) * t + sigma * W)
+  return(S)
+}
+# Function to calculate option price using Monte Carlo simulation
+monte_carlo_option_price <- function(S0, K, r, sigma, T, n, num_simulations) {
+  option_prices <- numeric(num_simulations)
+  for (i in 1:num_simulations) {
+    ST <- simulate_stock_prices(S0, r, sigma, T, n)[n + 1] # Final stock price
+    option_prices[i] <- pmax(ST - K, 0) # Payoff of the option
+  }
+  option_price <- mean(option_prices) * exp(-r * T) # Discounted expected payoff
+  return(option_price)
+}
+# Number of Monte Carlo simulations
+num_simulations <- 10000
+option_price <- monte_carlo_option_price(S0, K, r, sigma, T, n, num_simulations)
+cat("Option price:", option_price, "\n")