notebooks /2018-05-09-gibbs-sampling.ipynb

[waudinio27]

Hello Jessy! How are you? It is me Matthias. I found your blog and I think it is really great. I immediately liked it. The Gibbs part took so much my interest that I made a program in the same style inspired by you. The problem is that the predictions I added are not stable. Can you have a look and maybe see what can be the cause?

Thanks and best regards

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import yfinance as yf

# Fetch BTC-USD Data from yfinance
def fetch_yfinance_data(ticker, interval='1h'):
    # Download the data using yfinance (max available period)
    df = yf.download(ticker, interval=interval, period="max")
    
    # We're interested in 'Close' and 'Volume'
    df = df[['Close', 'Volume']].dropna()
    
    return df

# Create the joint distribution based on financial data
def setup_joint_distribution(data):
    # For simplicity, we model 'close' and 'volume' as Gaussian processes with covariance
    D = 2  # Two dimensions: close price and volume
    close = data['Close'].values
    volume = data['Volume'].values
    
    # Calculate means and covariance from data
    joint_mu = np.array([np.mean(close), np.mean(volume)])
    joint_cov = np.cov(close, volume)
    
    return joint_mu, joint_cov

# Set up conditional distribution like in the original code
def get_conditional_dist(joint_mu, joint_cov, var_index):
    a = joint_mu[var_index]
    b = joint_mu[1 - var_index]
    
    A = joint_cov[var_index, var_index]
    B = joint_cov[1 - var_index, 1 - var_index]
    C = joint_cov[var_index, 1 - var_index]
    
    B_inv = 1 / B

    def dist(g):
        mu = a + C * B_inv * (g - b)
        cov = A - B_inv * C * C
        return np.sqrt(cov) * np.random.randn(1) + mu
    
    return dist

# Gibbs Sampling for two financial variables (close and volume) and extend with forecast
def gibbs_sample_with_forecast(univariate_conditionals, sample_count=10000, forecast_steps=10):
    D = len(univariate_conditionals)
    samples = np.zeros((D, sample_count + forecast_steps))
    samples[:, 0] = [3, -3]  # Initialize with arbitrary values
    
    # Perform Gibbs sampling for historical data
    for i in range(1, sample_count):
        samples[:, i] = samples[:, i - 1]
        d = i % D
        samples[d, i] = univariate_conditionals[d](samples[1 - d, i - 1])
    
    # Forecast future values by extending the Gibbs sampling process
    for step in range(forecast_steps):
        i = sample_count + step
        samples[:, i] = samples[:, i - 1]
        d = step % D  # Alternate between close and volume
        samples[d, i] = univariate_conditionals[d](samples[1 - d, i - 1])
    
    return samples

# Visualization including forecasted future samples
def plot_samples_and_forecast(samples, true_distribution, forecasted_samples, forecast_steps):
    fig, ax = plt.subplots()
    ax.plot(*true_distribution, '.', alpha=0.1)
    ax.plot(*samples[:, :-forecast_steps], 'k')
    ax.plot(*samples[:, :-forecast_steps], '.r')
    
    # Plot forecasted samples
    future_x = np.arange(len(samples[0]) - forecast_steps, len(samples[0]))
    ax.plot(future_x, forecasted_samples[0, -forecast_steps:], 'b--', label='Forecast Close')
    ax.plot(future_x, forecasted_samples[1, -forecast_steps:], 'g--', label='Forecast Volume')
    ax.plot(future_x, forecasted_samples[0, -forecast_steps:], 'bo')
    ax.plot(future_x, forecasted_samples[1, -forecast_steps:], 'go')
    ax.legend()
    ax.axis('square')
    plt.show()

# Function to print forecast values for close and volume
def print_forecast_values(samples, forecast_steps):
    # Extract the forecasted "close" and "volume" values
    forecasted_close = samples[0, -forecast_steps:]
    forecasted_volume = samples[1, -forecast_steps:]
    
    # Print the forecasted values
    print("\nForecasted Close Prices:")
    for i, value in enumerate(forecasted_close, 1):
        print(f"Step {i}: {value}")
    
    print("\nForecasted Volume:")
    for i, value in enumerate(forecasted_volume, 1):
        print(f"Step {i}: {value}")

# Main process

def main():
    ticker = "BTC-USD"  # Use BTC-USD for Bitcoin data
    data = fetch_yfinance_data(ticker, interval='1h')  # Hourly data with max period

    # Set up joint distribution for close price and volume
    joint_mu, joint_cov = setup_joint_distribution(data)
    
    # Create conditional distributions for Gibbs sampling
    univariate_conditionals = [get_conditional_dist(joint_mu, joint_cov, d) for d in range(2)]
    
    # Perform Gibbs sampling and forecast next 10 steps
    sample_count = 100000
    forecast_steps = 5
    samples = gibbs_sample_with_forecast(univariate_conditionals, sample_count=sample_count, forecast_steps=forecast_steps)
    
    # Plot the sampled data and forecasted future values
    plot_samples_and_forecast(samples, samples, samples, forecast_steps)
    
    # Print forecasted close and volume values
    print_forecast_values(samples, forecast_steps)

if __name__ == "__main__":
    main()