{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":102335,"databundleVersionId":12518947,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Editting\n* I am trying statistical models—such as SARIMA—for feature engineering and predicting whether a transition occurs.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom statsmodels.tsa.arima.model import ARIMA\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nfrom statsmodels.tsa.stattools import adfuller, kpss\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\nfrom statsmodels.stats.diagnostic import acorr_ljungbox\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\nimport warnings\nwarnings.filterwarnings('ignore')\n\nRAW_DIR = Path(\"/kaggle/input/cmi-detect-behavior-with-sensor-data\")\ndf = pd.read_csv(RAW_DIR / \"train.csv\")\n\ndf_sample=df[df['sequence_id']==\"SEQ_000007\"]['acc_x']","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:11.812305Z","iopub.execute_input":"2025-07-05T16:43:11.812648Z","iopub.status.idle":"2025-07-05T16:43:45.931494Z","shell.execute_reply.started":"2025-07-05T16:43:11.812625Z","shell.execute_reply":"2025-07-05T16:43:45.930613Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Basic Statistics of the Data:\")\nprint(df_sample.describe())\nprint(f\"\\nNumber of Data Points: {len(df_sample)}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:45.932778Z","iopub.execute_input":"2025-07-05T16:43:45.933073Z","iopub.status.idle":"2025-07-05T16:43:45.941921Z","shell.execute_reply.started":"2025-07-05T16:43:45.933022Z","shell.execute_reply":"2025-07-05T16:43:45.941182Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 1. Initial time series visualization","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(2, 2, figsize=(15, 10))\nfig.suptitle('Initial Time Series Analysis', fontsize=16)\n\n# Original data plot\naxes[0, 0].plot(df_sample.index, df_sample.values, 'b-', linewidth=1)\naxes[0, 0].set_title('Original Time Series Data (acc_x)')\naxes[0, 0].set_xlabel('Time')\naxes[0, 0].set_ylabel('Acceleration (acc_x)')\naxes[0, 0].grid(True, alpha=0.3)\n\n# Histogram of values\naxes[0, 1].hist(df_sample.values, bins=30, alpha=0.7, edgecolor='black')\naxes[0, 1].set_title('Data Distribution')\naxes[0, 1].set_xlabel('Acceleration Values')\naxes[0, 1].set_ylabel('Frequency')\naxes[0, 1].grid(True, alpha=0.3)\n\n# Autocorrelation function (ACF)\nplot_acf(df_sample, ax=axes[1, 0], lags=40, title='Autocorrelation Function (ACF)')\naxes[1, 0].grid(True, alpha=0.3)\n\n# Partial autocorrelation function (PACF)\nplot_pacf(df_sample, ax=axes[1, 1], lags=27, title='Partial Autocorrelation Function (PACF)')\naxes[1, 1].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:45.943006Z","iopub.execute_input":"2025-07-05T16:43:45.943308Z","iopub.status.idle":"2025-07-05T16:43:47.075627Z","shell.execute_reply.started":"2025-07-05T16:43:45.943288Z","shell.execute_reply":"2025-07-05T16:43:47.074522Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 2. Stationarity test function","metadata":{}},{"cell_type":"code","source":"def check_stationarity(ts, title):\n    print(f\"\\n=== Stationarity Test for {title} ===\")\n    \n    adf_result = adfuller(ts, autolag='AIC')\n    print(\"ADF Test:\")\n    print(f\"  ADF Statistic: {adf_result[0]:.6f}\")\n    print(f\"  p-value: {adf_result[1]:.6f}\")\n    print(f\"  Critical Values: {adf_result[4]}\")\n    if adf_result[1] <= 0.05:\n        print(\"  Result: Stationary (p < 0.05)\")\n    else:\n        print(\"  Result: Non-stationary (p >= 0.05)\")\n    \n    kpss_result = kpss(ts, regression='c', nlags=\"auto\")\n    print(\"\\nKPSS Test:\")\n    print(f\"  KPSS Statistic: {kpss_result[0]:.6f}\")\n    print(f\"  p-value: {kpss_result[1]:.6f}\")\n    print(f\"  Critical Values: {kpss_result[3]}\")\n    if kpss_result[1] >= 0.05:\n        print(\"  Result: Stationary (p >= 0.05)\")\n    else:\n        print(\"  Result: Non-stationary (p < 0.05)\")\n\n# Run stationarity tests on the original data\ncheck_stationarity(df_sample, \"Original Data\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:47.077793Z","iopub.execute_input":"2025-07-05T16:43:47.078126Z","iopub.status.idle":"2025-07-05T16:43:47.096290Z","shell.execute_reply.started":"2025-07-05T16:43:47.078101Z","shell.execute_reply":"2025-07-05T16:43:47.095369Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 3. Estimate differencing order\n","metadata":{}},{"cell_type":"code","source":"def estimate_differencing_order(ts, max_d=2):\n    print(\"\\n=== Estimating Differencing Order ===\")\n    print(\"Note: SARIMA handles differencing internally; actual differencing is not applied here.\")\n    \n    current_ts = ts.copy()\n    recommended_d = 0\n    \n    for i in range(max_d + 1):\n        adf_stat, p_value = adfuller(current_ts, autolag='AIC')[:2]\n        print(f\"Differencing order {i}: ADF p-value = {p_value:.6f}\")\n        if p_value <= 0.05:\n            print(f\"Recommended differencing order: d = {i}\")\n            recommended_d = i\n            break\n        if i < max_d:\n            current_ts = current_ts.diff().dropna()\n            recommended_d = i + 1\n    \n    if recommended_d > max_d:\n        print(f\"Data is not stationary even after maximum differencing order: d = {max_d}\")\n        recommended_d = max_d\n    return recommended_d\n\nrecommended_d = estimate_differencing_order(df_sample)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:47.097315Z","iopub.execute_input":"2025-07-05T16:43:47.097608Z","iopub.status.idle":"2025-07-05T16:43:47.128726Z","shell.execute_reply.started":"2025-07-05T16:43:47.097584Z","shell.execute_reply":"2025-07-05T16:43:47.127874Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 4. Seasonality detection","metadata":{}},{"cell_type":"code","source":"def detect_seasonality(ts, max_lag=50):\n    print(\"\\n=== Detecting Seasonality ===\")\n    \n    fig, ax = plt.subplots(figsize=(12, 6))\n    plot_acf(ts, ax=ax, lags=max_lag, title='ACF for Seasonality Detection')\n    plt.show()\n\n    acf_values = [abs(ts.autocorr(lag=lag)) for lag in range(1, min(max_lag+1, len(ts)//4))]\n    lags = list(range(1, len(acf_values)+1))\n    significant_lags = [(lags[i], acf_values[i]) for i in range(len(acf_values)) if acf_values[i] > 0.1]\n    print(f\"Lags with significant autocorrelation: {significant_lags}\")\n\n    common_periods = [12, 24, 48, 168]\n    period_correlations = []\n    for period in common_periods:\n        if period <= len(ts)//4:\n            corr = abs(ts.autocorr(lag=period))\n            period_correlations.append((period, corr))\n            print(f\"Autocorrelation at period {period}: {corr:.4f}\")\n\n    max_corr_lag, max_corr_value = max(zip(lags, acf_values), key=lambda x: x[1])\n    if period_correlations:\n        best_common_period, best_common_corr = max(period_correlations, key=lambda x: x[1])\n    else:\n        best_common_period, best_common_corr = None, 0\n\n    if best_common_corr > 0.15 and best_common_corr > max_corr_value:\n        selected_period = best_common_period\n        print(f\"\\nSelected seasonal period: {selected_period} (common period, correlation: {best_common_corr:.4f})\")\n    elif max_corr_value > 0.1:\n        selected_period = max_corr_lag\n        print(f\"\\nSelected seasonal period: {selected_period} (highest autocorrelation lag, value: {max_corr_value:.4f})\")\n    else:\n        selected_period = 24 if len(ts) >= 24 else max(12, len(ts)//10)\n        print(f\"\\nNo clear seasonality detected; defaulting to: {selected_period}\")\n\n    # Optional FFT-based detection\n    try:\n        from scipy.fft import fft, fftfreq\n        fft_vals = np.abs(fft(ts.values))\n        fft_freqs = fftfreq(len(ts), 1)\n        pos = fft_freqs > 0\n        pos_fft = fft_vals[pos]\n        pos_freqs = fft_freqs[pos]\n        idx = np.argmax(pos_fft)\n        dominant_freq = pos_freqs[idx]\n        fft_period = int(1 / dominant_freq)\n        print(f\"FFT-based primary period: {fft_period}\")\n        if 6 <= fft_period <= len(ts)//4 and pos_fft[idx] > np.mean(pos_fft)*2:\n            print(\"Considering FFT result as well.\")\n    except ImportError:\n        print(\"Skipping FFT analysis.\")\n\n    return selected_period\n\nseasonal_period = detect_seasonality(df_sample)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:47.129819Z","iopub.execute_input":"2025-07-05T16:43:47.130150Z","iopub.status.idle":"2025-07-05T16:43:47.418093Z","shell.execute_reply.started":"2025-07-05T16:43:47.130122Z","shell.execute_reply":"2025-07-05T16:43:47.417101Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 5. Automatic SARIMA parameter selection based on AIC","metadata":{}},{"cell_type":"code","source":"def auto_sarima_selection(ts, seasonal_period, max_p=3, max_d=2, max_q=3, max_P=2, max_D=1, max_Q=2):\n    \"\"\"\n    Automatically select SARIMA parameters based on AIC.\n    Uses the detected seasonal period.\n    \"\"\"\n    print(\"\\n=== Automatic Selection of SARIMA Parameters ===\")\n    print(f\"Seasonal period to use: {seasonal_period}\")\n    \n    best_aic = float('inf')\n    best_params = None\n    best_seasonal_params = None\n    results = []\n    \n    # Try multiple candidate seasonal periods:\n    # Use both the detected period and common periods if data length allows\n    seasonal_candidates = [seasonal_period]\n    \n    common_periods = [12, 24, 168]  # monthly, daily, weekly\n    for period in common_periods:\n        if period <= len(ts) // 4 and period != seasonal_period:\n            seasonal_candidates.append(period)\n    \n    print(f\"Seasonal periods to try: {seasonal_candidates}\")\n    \n    for s in seasonal_candidates:\n        print(f\"\\nOptimizing with seasonal period {s}...\")\n        \n        # Iterate over parameter combinations\n        for p in range(max_p + 1):\n            for d in range(max_d + 1):\n                for q in range(max_q + 1):\n                    for P in range(max_P + 1):\n                        for D in range(max_D + 1):\n                            for Q in range(max_Q + 1):\n                                try:\n                                    model = SARIMAX(ts, \n                                                   order=(p, d, q),\n                                                   seasonal_order=(P, D, Q, s),\n                                                   enforce_stationarity=False,\n                                                   enforce_invertibility=False)\n                                    fitted_model = model.fit(disp=False)\n                                    \n                                    aic = fitted_model.aic\n                                    results.append({\n                                        'order': (p, d, q),\n                                        'seasonal_order': (P, D, Q, s),\n                                        'AIC': aic,\n                                        'BIC': fitted_model.bic\n                                    })\n                                    \n                                    if aic < best_aic:\n                                        best_aic = aic\n                                        best_params = (p, d, q)\n                                        best_seasonal_params = (P, D, Q, s)\n                                        \n                                except:\n                                    continue\n    \n    # Display the top 10 models by AIC\n    results_df = pd.DataFrame(results)\n    results_df = results_df.sort_values('AIC').head(10)\n    \n    print(\"\\nTop 10 models (sorted by AIC):\")\n    print(results_df.to_string(index=False))\n    \n    print(f\"\\nBest parameters:\")\n    print(f\"  ARIMA order: {best_params}\")\n    print(f\"  Seasonal order: {best_seasonal_params}\")\n    print(f\"  AIC: {best_aic:.4f}\")\n    \n    return best_params, best_seasonal_params, best_aic\n\n# Automatic parameter selection (using detected seasonal period)\nbest_order, best_seasonal_order, best_aic = auto_sarima_selection(\n    df_sample, seasonal_period, max_p=2, max_d=1, max_q=2, max_P=1, max_D=1, max_Q=1\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:43:47.419104Z","iopub.execute_input":"2025-07-05T16:43:47.419435Z","iopub.status.idle":"2025-07-05T16:44:17.972416Z","shell.execute_reply.started":"2025-07-05T16:43:47.419408Z","shell.execute_reply":"2025-07-05T16:44:17.971549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Fitting the Optimal SARIMA Model ===\nprint(\"\\n=== Fitting the Optimal SARIMA Model ===\")\n\n# Create and fit the model\nmodel = SARIMAX(df_sample, \n                order=best_order,\n                seasonal_order=best_seasonal_order,\n                enforce_stationarity=False,\n                enforce_invertibility=False)\n\nfitted_model = model.fit(disp=False)\n\n# Display model summary\nprint(fitted_model.summary())\n\n# === SARIMA Model Diagnostics ===\nfig, axes = plt.subplots(2, 2, figsize=(15, 10))\nfig.suptitle('SARIMA Model Diagnostics', fontsize=16)\n\n# Residual Partial Autocorrelation (PACF)\nresiduals = fitted_model.resid\nplot_pacf(residuals, ax=axes[0, 0], lags=27, title='Residual Partial Autocorrelation (PACF)')\naxes[0, 0].grid(True, alpha=0.3)\n\n# Residual Autocorrelation (ACF)\nplot_acf(residuals, ax=axes[0, 1], lags=40, title='Residual Autocorrelation (ACF)')\naxes[0, 1].grid(True, alpha=0.3)\n\n# Residual Time Series Plot\naxes[1, 0].plot(residuals.index, residuals.values, 'r-', linewidth=1)\naxes[1, 0].set_title('Residual Time Series Plot')\naxes[1, 0].set_xlabel('Time')\naxes[1, 0].set_ylabel('Residuals')\naxes[1, 0].grid(True, alpha=0.3)\n\n# Residual Distribution\naxes[1, 1].hist(residuals.values, bins=30, alpha=0.7, edgecolor='black')\naxes[1, 1].set_title('Residual Distribution')\naxes[1, 1].set_xlabel('Residuals')\naxes[1, 1].set_ylabel('Frequency')\naxes[1, 1].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:44:17.973293Z","iopub.execute_input":"2025-07-05T16:44:17.973511Z","iopub.status.idle":"2025-07-05T16:44:19.182486Z","shell.execute_reply.started":"2025-07-05T16:44:17.973494Z","shell.execute_reply":"2025-07-05T16:44:19.181543Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === SARIMA Model Forecast Results ===\nfig, axes = plt.subplots(2, 1, figsize=(15, 10))\nfig.suptitle('SARIMA Model Forecast Results', fontsize=16)\n\n# Fitted values and confidence intervals\nfitted_values = fitted_model.fittedvalues\nconf_int = fitted_model.get_forecast(steps=len(df_sample)).conf_int()\n\n\n# Detailed view for the last n_zoom points\nn_zoom = min(57, len(df_sample))\nzoom_idx = df_sample.index[-n_zoom:]\naxes[0].plot(zoom_idx, df_sample.loc[zoom_idx], 'b-', linewidth=2, label='Actual Data', alpha=0.7)\naxes[0].plot(zoom_idx, fitted_values.loc[zoom_idx], 'r-', linewidth=2, label='SARIMA Fit', alpha=0.8)\naxes[0].fill_between(zoom_idx,\n                     fitted_values.loc[zoom_idx] - 1.96*np.std(residuals),\n                     fitted_values.loc[zoom_idx] + 1.96*np.std(residuals),\n                     alpha=0.2, label='95% Confidence Interval')\naxes[0].set_title(f'Detailed View (Last {n_zoom} Data Points')\naxes[0].set_xlabel('Time')\naxes[0].set_ylabel('Acceleration (acc_x)')\naxes[0].legend()\naxes[0].grid(True, alpha=0.3)\n\n# Detailed view for the last n_zoom points\nn_zoom = min(30, len(df_sample))\nzoom_idx = df_sample.index[-n_zoom:]\naxes[1].plot(zoom_idx, df_sample.loc[zoom_idx], 'b-', linewidth=2, label='Actual Data', alpha=0.7)\naxes[1].plot(zoom_idx, fitted_values.loc[zoom_idx], 'r-', linewidth=2, label='SARIMA Fit', alpha=0.8)\naxes[1].fill_between(zoom_idx,\n                     fitted_values.loc[zoom_idx] - 1.96*np.std(residuals),\n                     fitted_values.loc[zoom_idx] + 1.96*np.std(residuals),\n                     alpha=0.2, label='95% Confidence Interval')\naxes[1].set_title(f'Detailed View (Last {n_zoom} Data Points')\naxes[1].set_xlabel('Time')\naxes[1].set_ylabel('Acceleration (acc_x)')\naxes[1].legend()\naxes[1].grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:44:19.183534Z","iopub.execute_input":"2025-07-05T16:44:19.183797Z","iopub.status.idle":"2025-07-05T16:44:19.861561Z","shell.execute_reply.started":"2025-07-05T16:44:19.183777Z","shell.execute_reply":"2025-07-05T16:44:19.860636Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Model Performance Evaluation ===\nprint(\"\\n=== Model Performance Evaluation ===\")\n\n# Basic evaluation metrics\nmse = mean_squared_error(df_sample, fitted_values)\nrmse = np.sqrt(mse)\nmae = mean_absolute_error(df_sample, fitted_values)\nmape = np.mean(np.abs((df_sample - fitted_values) / df_sample)) * 100\n\nprint(f\"Mean Squared Error (MSE): {mse:.6f}\")\nprint(f\"Root Mean Squared Error (RMSE): {rmse:.6f}\")\nprint(f\"Mean Absolute Error (MAE): {mae:.6f}\")\nprint(f\"Mean Absolute Percentage Error (MAPE): {mape:.2f}%\")\n\n# Residual Normality Test (Shapiro-Wilk)\nfrom scipy import stats\nshapiro_stat, shapiro_p = stats.shapiro(residuals)\nprint(f\"\\nResidual Normality Test (Shapiro-Wilk):\")\nprint(f\"  Statistic: {shapiro_stat:.6f}\")\nprint(f\"  p-value: {shapiro_p:.6f}\")\n\n# Residual Independence Test (Ljung-Box)\nljung_box_result = acorr_ljungbox(residuals, lags=10, return_df=True)\nprint(f\"\\nResidual Independence Test (Ljung-Box):\")\nprint(ljung_box_result)\n\n# === Future Forecast ===\n\n# Forecast for the next 30 periods\nn_forecast = 30\nforecast = fitted_model.get_forecast(steps=n_forecast)\nforecast_values = forecast.predicted_mean\nforecast_conf_int = forecast.conf_int()\n\n# Create index for forecast period\nlast_step = df_sample.index[-1]  # これは整数\nforecast_index = pd.Index(\n    np.arange(last_step + 1, last_step + 1 + n_forecast),\n    name=df_sample.index.name\n)\n# Plot future forecast\nfig, ax = plt.subplots(figsize=(15, 8))\n\n# Last n_history points of historical data\nn_history = min(50, len(df_sample))\nhistory_idx = df_sample.index[-n_history:]\n\nax.plot(history_idx, df_sample.loc[history_idx], 'b-', linewidth=2, label='Actual Data', alpha=0.7)\nax.plot(forecast_index, forecast_values, 'r-', linewidth=2, label='Forecast', alpha=0.8)\nax.fill_between(forecast_index,\n                forecast_conf_int.iloc[:, 0],\n                forecast_conf_int.iloc[:, 1],\n                alpha=0.2, label='95% Confidence Interval')\n\nax.axvline(x=last_step, linestyle='--', alpha=0.7, label='Forecast Start')\nax.set_title('Future Forecast by SARIMA Model')\nax.set_xlabel('Time')\nax.set_ylabel('Acceleration (acc_x)')\nax.legend()\nax.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.show()\n\nprint(f\"\\nForecast Statistics:\")\nprint(f\"  Forecast Mean: {forecast_values.mean():.6f}\")\nprint(f\"  Forecast Std Dev: {forecast_values.std():.6f}\")\nprint(f\"  Forecast Min: {forecast_values.min():.6f}\")\nprint(f\"  Forecast Max: {forecast_values.max():.6f}\")\n\nprint(\"\\n=== Analysis Complete ===\")\nprint(\"SARIMA model application is complete.\")\nprint(\"Parameter selection was carried out with the following steps:\")\nprint(\"1. Stationarity tests (ADF, KPSS tests)\")\nprint(\"2. Estimation of differencing orders (actual differencing not applied since SARIMA processes internally)\")\nprint(\"3. Seasonal detection (identifying seasonal periods using multiple methods)\")\nprint(\"4. Automatic parameter selection based on AIC\")\nprint(\"5. Model fitting and diagnostics\")\n\nprint(\"\\nKey Improvements:\")\nprint(\"- Leveraged SARIMA model's differencing capabilities\")\nprint(\"- Analysis preserving original data information\")\nprint(\"- Comparison of multiple seasonal period candidates\")\nprint(\"- Efficient search guided by recommended differencing orders\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-05T16:44:19.863467Z","iopub.execute_input":"2025-07-05T16:44:19.863715Z","iopub.status.idle":"2025-07-05T16:44:20.425081Z","shell.execute_reply.started":"2025-07-05T16:44:19.863695Z","shell.execute_reply":"2025-07-05T16:44:20.424228Z"}},"outputs":[],"execution_count":null}]}