moving_average_python_threaded.py

import numpy as np
import pandas as pd
from numba import jit
import time
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
import multiprocessing as mp
from functools import partial
import os


def calculate_basic_features_chunk(args):
    """Calculate basic features for a chunk of data"""
    opens, highs, lows, closes, volumes, start_idx = args
    n = len(closes)
    chunk_size = len(closes)
    features = np.zeros((chunk_size, 50))  # Start with 50 basic features
    
    for i in range(chunk_size):
        global_i = start_idx + i
        # Price-based features
        features[i, 0] = closes[i]  # Close price
        features[i, 1] = opens[i]   # Open price
        features[i, 2] = highs[i]   # High price
        features[i, 3] = lows[i]    # Low price
        
        # Basic ratios and differences
        features[i, 4] = (closes[i] - opens[i]) / opens[i] if opens[i] != 0 else 0  # Return
        features[i, 5] = (highs[i] - lows[i]) / opens[i] if opens[i] != 0 else 0   # True range
        features[i, 6] = (closes[i] - lows[i]) / (highs[i] - lows[i]) if (highs[i] - lows[i]) != 0 else 0.5  # Stochastic
        features[i, 7] = volumes[i] if i < len(volumes) else 0  # Volume
        
        # Simple moving statistics (if enough history)
        if global_i >= 5:
            recent_close = closes[max(0, i-5):i+1]  # Note: using i instead of global_i for indexing within chunk
            if len(recent_close) > 0:
                features[i, 8] = np.mean(recent_close)  # 5-period SMA
                features[i, 9] = np.std(recent_close)   # 5-period volatility
                if np.max(recent_close) != np.min(recent_close):
                    features[i, 10] = (closes[i] - np.min(recent_close)) / (np.max(recent_close) - np.min(recent_close))  # 5-period percentile
                else:
                    features[i, 10] = 0.5
        
        if global_i >= 10:
            recent_close_10 = closes[max(0, i-10):i+1]
            if len(recent_close_10) > 0:
                features[i, 11] = np.mean(recent_close_10)  # 10-period SMA
                features[i, 12] = np.std(recent_close_10)   # 10-period volatility
                if np.max(recent_close_10) != np.min(recent_close_10):
                    features[i, 13] = (closes[i] - np.min(recent_close_10)) / (np.max(recent_close_10) - np.min(recent_close_10))  # 10-period percentile
                else:
                    features[i, 13] = 0.5
        
        if global_i >= 20:
            recent_close_20 = closes[max(0, i-20):i+1]
            if len(recent_close_20) > 0:
                features[i, 14] = np.mean(recent_close_20)  # 20-period SMA
                features[i, 15] = np.std(recent_close_20)   # 20-period volatility
                if np.max(recent_close_20) != np.min(recent_close_20):
                    features[i, 16] = (closes[i] - np.min(recent_close_20)) / (np.max(recent_close_20) - np.min(recent_close_20))  # 20-period percentile
                else:
                    features[i, 16] = 0.5
        
        # Momentum indicators
        if i >= 1 and global_i >= 1:
            features[i, 17] = (closes[i] - closes[i-1]) / closes[i-1] if closes[i-1] != 0 else 0  # 1-period return
        if i >= 3 and global_i >= 3:
            features[i, 18] = (closes[i] - closes[i-3]) / closes[i-3] if closes[i-3] != 0 else 0  # 3-period return
        if i >= 5 and global_i >= 5:
            features[i, 19] = (closes[i] - closes[i-5]) / closes[i-5] if closes[i-5] != 0 else 0  # 5-period return
        
        # Price position relative to recent highs/lows
        if global_i >= 10:
            recent_high_10 = highs[max(0, i-10):i+1]
            recent_low_10 = lows[max(0, i-10):i+1]
            if len(recent_high_10) > 0 and len(recent_low_10) > 0:
                recent_max = np.max(recent_high_10)
                recent_min = np.min(recent_low_10)
                if recent_max != recent_min:
                    features[i, 20] = (closes[i] - recent_min) / (recent_max - recent_min)  # Position in 10-day range
                else:
                    features[i, 20] = 0.5
                if recent_max != recent_min:
                    features[i, 21] = (highs[i] - recent_min) / (recent_max - recent_min)  # High position in 10-day range
                else:
                    features[i, 21] = 0.5
                if recent_max != recent_min:
                    features[i, 22] = (lows[i] - recent_min) / (recent_max - recent_min)  # Low position in 10-day range
                else:
                    features[i, 22] = 0.5
        
        # Volatility measures
        if global_i >= 10:
            returns = np.diff(closes[max(0, i-10):i+1])
            if len(returns) > 1:
                features[i, 23] = np.std(returns) if len(returns) > 1 else 0  # 10-period return volatility
                features[i, 24] = np.mean(returns)  # 10-period average return
        
        # Range-based features
        features[i, 25] = (highs[i] - lows[i]) / opens[i] if opens[i] != 0 else 0  # Daily range
        features[i, 26] = (highs[i] - closes[i]) / opens[i] if opens[i] != 0 else 0  # Upper shadow
        features[i, 27] = (closes[i] - lows[i]) / opens[i] if opens[i] != 0 else 0  # Lower shadow
        features[i, 28] = abs(opens[i] - closes[i]) / opens[i] if opens[i] != 0 else 0  # Body size
        
        # Log returns
        if i >= 1 and global_i >= 1 and closes[i-1] != 0:
            features[i, 29] = np.log(closes[i] / closes[i-1])
        else:
            features[i, 29] = 0
        
        # Higher moment features (manual calculations)
        if global_i >= 20:
            recent_returns = np.diff(closes[max(0, i-20):i+1])
            if len(recent_returns) > 2:
                mean_ret = np.mean(recent_returns)
                std_ret = np.std(recent_returns)
                if std_ret != 0:
                    # Manual skewness calculation
                    n_rets = len(recent_returns)
                    skew_num = np.sum(((recent_returns - mean_ret) / std_ret) ** 3)
                    features[i, 30] = skew_num * n_rets / ((n_rets - 1) * (n_rets - 2))  # Adjusted skewness
                    # Manual kurtosis calculation
                    kurt = np.mean(((recent_returns - mean_ret) / std_ret) ** 4)
                    features[i, 31] = kurt - 3  # Excess kurtosis
                else:
                    features[i, 30] = 0
                    features[i, 31] = 0
            else:
                features[i, 30] = 0
                features[i, 31] = 0
        
        # Volume features (if available)
        if i < len(volumes) and global_i >= 5:
            recent_volumes = volumes[max(0, i-5):i+1]
            features[i, 32] = volumes[i] / np.mean(recent_volumes) if np.mean(recent_volumes) != 0 else 1  # Volume ratio to recent avg
            features[i, 33] = volumes[i]  # Absolute volume
        
        # RSI-like indicator
        if global_i >= 14:
            total_gain = 0.0
            total_loss = 0.0
            count = 0
            for j in range(max(0, i-13), i+1):
                if j > 0 and global_i-13+j >= 1:  # Adjust for global indexing
                    change = closes[j] - closes[j-1]
                    if change > 0:
                        total_gain += change
                    else:
                        total_loss += abs(change)
                    count += 1
            
            avg_gain = total_gain / count if count > 0 else 0
            avg_loss = total_loss / count if count > 0 else 0
            
            if avg_loss != 0:
                features[i, 34] = 100 - (100 / (1 + avg_gain / avg_loss))  # RSI approximation
            else:
                features[i, 34] = 100
        
        # Bollinger Bands components
        if global_i >= 20:
            recent_closes_bb = closes[max(0, i-20):i+1]
            if len(recent_closes_bb) > 0:
                bb_mean = np.mean(recent_closes_bb)
                bb_std = np.std(recent_closes_bb)
                if bb_std != 0:
                    features[i, 35] = (closes[i] - bb_mean) / bb_std  # Bollinger Band position
                else:
                    features[i, 35] = 0
        
        # Directional features
        features[i, 36] = 1 if closes[i] > opens[i] else 0  # Bullish/Bearish
        if i > 0 and global_i > 0:
            features[i, 37] = 1 if highs[i] > highs[i-1] else 0  # New high
            features[i, 38] = 1 if lows[i] < lows[i-1] else 0  # New low
        
        # Gap features
        if i > 0 and global_i > 0:
            features[i, 39] = (opens[i] - closes[i-1]) / closes[i-1] if closes[i-1] != 0 else 0  # Gap up/down
            features[i, 40] = 1 if opens[i] > closes[i-1] else 0  # Gap up indicator
            features[i, 41] = 1 if opens[i] < closes[i-1] else 0  # Gap down indicator
        
        # Price level features
        features[i, 42] = np.log(closes[i]) if closes[i] > 0 else 0  # Log price
        features[i, 43] = closes[i] - opens[i]  # Body difference
        features[i, 44] = (highs[i] - max(opens[i], closes[i])) / opens[i] if opens[i] != 0 else 0  # Upper shadow normalized
        features[i, 45] = (min(opens[i], closes[i]) - lows[i]) / opens[i] if opens[i] != 0 else 0  # Lower shadow normalized
        
        # Acceleration features
        if i >= 2 and global_i >= 2:
            if closes[i-2] != 0 and closes[i-1] != 0:
                prev_return = (closes[i-1] - closes[i-2]) / closes[i-2] if closes[i-2] != 0 else 0
                curr_return = (closes[i] - closes[i-1]) / closes[i-1] if closes[i-1] != 0 else 0
                features[i, 46] = curr_return - prev_return  # Return acceleration
        
        # Trend strength (manual linear regression)
        if global_i >= 10:
            recent_closes = closes[max(0, i-10):i+1]
            n_regr = len(recent_closes)
            if n_regr > 1:
                x = np.arange(n_regr)
                y = recent_closes

                sum_x = np.sum(x)
                sum_y = np.sum(y)
                sum_xy = np.sum(x * y)
                sum_x2 = np.sum(x * x)

                denominator = n_regr * sum_x2 - sum_x * sum_x
                if denominator != 0:
                    slope = (n_regr * sum_xy - sum_x * sum_y) / denominator
                    features[i, 47] = slope  # Trend direction and strength
                else:
                    features[i, 47] = 0
            else:
                features[i, 47] = 0
        
        # Support/resistance proxies
        if global_i >= 10:
            features[i, 48] = np.std(closes[max(0, i-10):i+1])  # Volatility as resistance strength proxy
            features[i, 49] = np.mean(np.abs(np.diff(closes[max(0, i-10):i+1]))) if len(closes[max(0, i-10):i+1]) > 1 else 0  # Mean absolute change
    
    return features


def calculate_additional_features_chunk(args):
    """Calculate additional technical analysis features for a chunk"""
    opens, highs, lows, closes, volumes, start_idx = args
    n = len(closes)
    additional_features = np.zeros((n, 51))  # Additional 51 features
    
    for i in range(n):
        global_i = start_idx + i
        
        if global_i >= 14:
            # RSI approximation (14-period)
            gains = 0.0
            losses = 0.0
            for j in range(max(1, i-13), i+1):  # Adjust for chunk indexing
                if j > 0 and global_i-13+j-1 >= 0:  # Ensure we don't go out of bounds
                    change = closes[j] - closes[j-1]
                    if change > 0:
                        gains += change
                    else:
                        losses += abs(change)
            
            avg_gain = gains / min(14, i+1) if i+1 > 0 else 0  # Adjust for actual available data
            avg_loss = losses / min(14, i+1) if i+1 > 0 else 0
            
            if avg_loss != 0:
                rs = avg_gain / avg_loss
                additional_features[i, 0] = 100 - (100 / (1 + rs))  # RSI
            else:
                additional_features[i, 0] = 50  # Neutral value when no losses
        
        if global_i >= 12:
            # Simple MACD approximation (12, 26 period)
            period_short = min(12, i+1)
            period_long = min(26, i+1)
            if period_short > 0:
                ema_short = np.mean(closes[max(0, i-period_short+1):i+1])  # Simplified EMA
            else:
                ema_short = 0
            if period_long > 0:
                ema_long = np.mean(closes[max(0, i-period_long+1):i+1])   # Simplified EMA
            else:
                ema_long = 0
            additional_features[i, 1] = ema_short - ema_long  # MACD line approximation
        
        # Volatility (ATR approximation)
        if global_i >= 14:
            tr_list = []
            for j in range(max(1, i-14), i+1):
                if j > 0 and global_i-14+j-1 >= 0:  # Ensure we don't go out of bounds
                    h_l = highs[j] - lows[j]
                    h_pc = abs(highs[j] - closes[j-1])
                    l_pc = abs(lows[j] - closes[j-1])
                    true_range = max(h_l, h_pc, l_pc)
                    tr_list.append(true_range)
            if tr_list:
                additional_features[i, 11] = np.mean(tr_list)  # ATR approximation
        
        # Volume indicators
        if i >= 5 and len(volumes) > i:
            additional_features[i, 8] = volumes[i]  # Volume
            additional_features[i, 9] = volumes[i] / np.mean(volumes[max(0, i-5):i+1]) if np.mean(volumes[max(0, i-5):i+1]) != 0 else 1  # Volume ratio
    
    return additional_features


def split_data_for_parallel_processing(opens, highs, lows, closes, volumes, num_chunks=None):
    """Split data into chunks for parallel processing"""
    if num_chunks is None:
        num_chunks = min(os.cpu_count(), 8)  # Limit to 8 chunks to avoid overhead
    
    chunk_size = len(closes) // num_chunks
    chunks = []
    
    for i in range(num_chunks):
        start_idx = i * chunk_size
        end_idx = (i + 1) * chunk_size if i < num_chunks - 1 else len(closes)
        
        chunk_opens = opens[start_idx:end_idx]
        chunk_highs = highs[start_idx:end_idx]
        chunk_lows = lows[start_idx:end_idx]
        chunk_closes = closes[start_idx:end_idx]
        chunk_volumes = volumes[start_idx:end_idx]
        
        chunks.append((chunk_opens, chunk_highs, chunk_lows, chunk_closes, chunk_volumes, start_idx))
    
    return chunks


def calculate_quantitative_features_parallel(opens, highs, lows, closes, volumes):
    """Calculate quantitative features using parallel processing"""
    print(f"Using {os.cpu_count()} CPU cores for parallel processing")
    
    # Split data into chunks
    chunks = split_data_for_parallel_processing(opens, highs, lows, closes, volumes)
    
    # Process chunks in parallel
    with ProcessPoolExecutor() as executor:
        basic_feature_chunks = list(executor.map(calculate_basic_features_chunk, chunks))
    
    # Process additional features in parallel
    chunks = split_data_for_parallel_processing(opens, highs, lows, closes, volumes)
    with ProcessPoolExecutor() as executor:
        additional_feature_chunks = list(executor.map(calculate_additional_features_chunk, chunks))
    
    # Combine results
    all_basic_features = np.vstack(basic_feature_chunks)
    all_additional_features = np.vstack(additional_feature_chunks)
    
    # Combine all features
    all_features = np.concatenate([all_basic_features, all_additional_features], axis=1)
    
    return all_features


def calculate_moving_average_chunk(args):
    """Calculate moving average for a chunk of data"""
    prices, period, start_idx, end_idx = args
    ma_values = []
    
    # Calculate MA for the specified range
    for i in range(start_idx, min(end_idx, len(prices) - period + 1)):
        window = prices[i:i + period]
        if len(window) == period:
            ma_values.append(np.mean(window))
    
    return ma_values, start_idx


def calculate_multiple_moving_averages_parallel(closes, periods, num_workers=None):
    """Calculate multiple moving averages in parallel"""
    if num_workers is None:
        num_workers = min(os.cpu_count(), 4)  # Use fewer workers for MA calculation to avoid overhead
    
    results = {}
    
    for period in periods:
        print(f"Calculating MA_{period} in parallel...")
        
        # Split the work for this period
        chunk_size = max(1, (len(closes) - period + 1) // num_workers)
        chunks = []
        
        for i in range(num_workers):
            start_idx = i * chunk_size
            end_idx = (i + 1) * chunk_size if i < num_workers - 1 else len(closes) - period + 1
            chunks.append((closes, period, start_idx, end_idx))
        
        # Process chunks in parallel
        with ProcessPoolExecutor(max_workers=num_workers) as executor:
            chunk_results = list(executor.map(calculate_moving_average_chunk, chunks))
        
        # Combine results
        all_ma_values = []
        for ma_vals, orig_start_idx in chunk_results:
            all_ma_values.extend(ma_vals)
        
        results[f'MA_{period}'] = np.array(all_ma_values)
    
    return results


def main():
    start_time = time.time()
    
    print("Reading CSV file with pandas...")
    df = pd.read_csv("USDJPY2.csv")
    print(f"Loaded {len(df)} records.")
    
    # Extract price arrays
    opens = df['Open'].values.astype(np.float64)
    highs = df['High'].values.astype(np.float64)
    lows = df['Low'].values.astype(np.float64)
    closes = df['Close'].values.astype(np.float64)
    volumes = df['volume'].values.astype(np.float64) if 'volume' in df.columns else np.ones(len(df)) * 1000
    
    print("Calculating 100+ quantitative features for each row using parallel processing...")
    features = calculate_quantitative_features_parallel(opens, highs, lows, closes, volumes)
    print(f"Calculated {features.shape[1]} quantitative features for {features.shape[0]} rows.")
    
    # Calculate moving averages for periods 200-220 using parallel processing
    print("Calculating moving averages for periods 200-220 using parallel processing...")
    ma_periods = list(range(200, 221))  # 200 to 220 inclusive
    
    all_mas = calculate_multiple_moving_averages_parallel(closes, ma_periods)
    
    for period, ma_values in all_mas.items():
        print(f"Calculated {len(ma_values)} {period} values.")
    
    end_time = time.time()
    duration = (end_time - start_time) * 1000  # Convert to milliseconds
    print(f"Total execution time: {duration:.2f} ms")
    print(f"Features shape: {features.shape}")


if __name__ == "__main__":
    main()