<jupyter_start><jupyter_text>Licensed under the MIT License.Copyright (c) 2021-2031. All rights reserved. Kats Outliers Detection* Kats General  * `TimeSeriesData` params and methods: https://facebookresearch.github.io/Kats/api/kats.consts.htmlkats.consts.TimeSeriesData* Kats Detection  * Kats detection official tutorial: https://github.com/facebookresearch/Kats/blob/main/tutorials/kats_202_detection.ipynb    * It describes Kats' outlier detector's algorithms    * But Kats' multivariate anomaly detection only output strange errors to me, even using the same tutorial code, see this ticket: https://github.com/facebookresearch/Kats/issues/194* Other Kats Outlier Detectors  * https://facebookresearch.github.io/Kats/api/kats.detectors.prophet_detector.html    * Kats v0.1 requires prophet version to be "0.7" exactly, other will get errors, but my laptop could only install higher version prophet...  * https://facebookresearch.github.io/Kats/api/kats.detectors.hourly_ratio_detection.html    * It requires the time series to be hour-level granularity<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import kpss
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.stattools import durbin_watson
from statsmodels.tsa.api import VAR
from statsmodels.tsa.vector_ar.vecm import VECM

from kats.consts import TimeSeriesData
from kats.detectors.outlier import OutlierDetector

import warnings
warnings.filterwarnings("ignore")
ts_df = pd.read_pickle('../../crystal_ball/data_collector/structured_data/sales_ts.pkl')
print(ts_df.shape)

ts_df_train = ts_df.iloc[ts_df.index < '2015-03-01']
print(ts_df_train.shape)

ts_df.head()
def plot_ts(ts, title):
    plt.figure(figsize=(20,3))
    for col in ts.columns:
        fig = plt.plot(ts[col], label=col)
    plt.title(title)
    plt.legend(loc='best')
    plt.tight_layout()
    plt.show()
    

def plot_ts_outliers(ts, title, outliers, decomp='additive'):
    outliers_x = [str(outlier).split()[0] for outlier in outliers[0]]
    outliers_y = ts.iloc[ts.index.isin(outliers_x)]
    
    plt.figure(figsize=(20,10))
    plt.subplot(411)
    fig = plt.plot(ts, label='original ts', color='blue')
    plt.scatter(outliers_x, outliers_y, c='red', marker='*')
    plt.legend(loc='best')
    
    plt.subplot(412)
    decomposition = seasonal_decompose(ts, model=decomp)
    residual = decomposition.resid
    fig = plt.plot(residual, label='residuals', color='purple')
    outliers_y_res = residual.iloc[residual.index.isin(outliers_x)]
    plt.scatter(outliers_x, outliers_y_res, c='red', marker='*')
    plt.legend(loc='best')
    
    plt.title(title)
    plt.tight_layout()
    plt.show()
plot_ts(ts_df_train, title='Univariate training ts plot')
# Covert to Kats required TimeSeriesData input

kats_ts_all = TimeSeriesData(ts_df_train.reset_index().rename(index=str, columns={'Date': 'time'}))
print(len(kats_ts_all))<jupyter_output>789<jupyter_text>Univariate OutlierDetector* Kats' outlier detector: https://facebookresearch.github.io/Kats/api/kats.detectors.outlier.html<jupyter_code># detect & plot outliers

ts_outlierDetection = OutlierDetector(kats_ts_all, 'multiplication', iqr_mult=5)
ts_outlierDetection.detector()

plot_ts_outliers(ts_df_train, title='Outliers in all ts train', outliers=ts_outlierDetection.outliers, decomp='multipllicative')
# remove and plot outliers

ts_outlierDetection_outliers_removed = ts_outlierDetection.remover(interpolate = False) # No interpolation
ts_outlierDetection_interpolated = ts_outlierDetection.remover(interpolate = True) # With linear interpolation

ts_outlierDetection_outliers_removed
fig, ax = plt.subplots(figsize=(25,8), nrows=1, ncols=2)
ts_outlierDetection_outliers_removed.to_dataframe().plot(x='time',y = 'y_0', ax = ax[0])
ax[0].set_title("Outliers Removed : No interpolation")
ts_outlierDetection_interpolated.to_dataframe().plot(x = 'time',y = 'y_0', ax = ax[1])
ax[1].set_title("Outliers Removed : With interpolation")
plt.show()
sub_original_df = ts_df_train.iloc[(ts_df_train.index>='2013-12-22') & (ts_df_train.index<='2014-01-02')]
sub_df_removed = ts_outlierDetection_outliers_removed.to_dataframe()
sub_df_removed = sub_df_removed.loc[(sub_df_removed['time']>='2013-12-22') & (sub_df_removed['time']<='2014-01-02')]

sub_df_interpolated = ts_outlierDetection_interpolated.to_dataframe()
sub_df_interpolated = sub_df_interpolated.loc[(sub_df_interpolated['time']>='2013-12-22') & (sub_df_interpolated['time']<='2014-01-02')]

fig, ax = plt.subplots(figsize=(25,8), nrows=1, ncols=2)
sub_original_df.reset_index().plot(x='Date', y='Daily_Sales', ax=ax[0], color='orange', marker='o', label='original ts')
sub_df_removed.plot(x='time', y='y_0', ax= ax[0], color='green', label='outlier removed ts')
ax[0].set_title("Outliers Removed Subset: No interpolation")

sub_original_df.reset_index().plot(x='Date', y='Daily_Sales', ax=ax[1], color='orange', marker='o', label='original ts')
sub_df_interpolated.plot(x = 'time',y = 'y_0', ax= ax[1], color='green', label='outlier interpolated ts')
ax[1].set_title("Outliers Removed Subset: With interpolation")
plt.show()<jupyter_output><empty_output><jupyter_text>Multivariate Anomaly Detection* References  * VAR for anomaly detection: https://www.analyticsvidhya.com/blog/2021/08/multivariate-time-series-anomaly-detection-using-var-model/  * More about VAR: https://www.machinelearningplus.com/time-series/vector-autoregression-examples-python/* More multivariate time series models: https://www.statsmodels.org/dev/api.htmlmultivariate-time-series-models<jupyter_code>mul_ts_df = pd.read_pickle('../../crystal_ball/data_collector/structured_data/multivar_ts.pkl')

print(mul_ts_df.shape)
mul_ts_df.head()
occupancy = mul_ts_df[['Occupancy']]
mul_ts_df.drop('Occupancy', inplace=True, axis=1)
print(mul_ts_df.shape)<jupyter_output>(8143, 5)<jupyter_text>Convert Data to Stationary<jupyter_code>def test_stationarity_multi_ts(multi_ts_df):
    results_dct = {}
    
    for col in multi_ts_df.columns:
        timeseries = multi_ts_df[col]
        
        adf_result, kpss_result = None, None
        results_dct[col] = {'Differencing Stationary': None, 'Trending Stationary': None}
    
        # Perform Augmented Dickey-Fuller test:
        adftest = adfuller(timeseries, autolag='AIC')
        adf_output = pd.Series(adftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
        adf_test_stats = adf_output['Test Statistic']
        for key,value in adftest[4].items():
            adf_output[f'Critical Value {key}'] = value
            
        if abs(adf_test_stats) >= abs(adf_output[f'Critical Value 1%']):
            adf_result = '99%'
        elif abs(adf_test_stats) >= abs(adf_output[f'Critical Value 5%']) and abs(adf_test_stats) < abs(adf_output[f'Critical Value 1%']):
            adf_result = '95%'
        elif abs(adf_test_stats) >= abs(adf_output[f'Critical Value 10%']) and abs(adf_test_stats) < abs(adf_output[f'Critical Value 5%']):
            adf_result = '90%'

            
        # Perform KPSS
        kpsstest = kpss(timeseries, regression='c')
        kpss_output = pd.Series(kpsstest[0:3], index=['Test Statistic','p-value','Lags Used'])
        kpss_test_stats = kpss_output['Test Statistic']
        for key,value in kpsstest[3].items():
            kpss_output[f'Critical Value {key}'] = value

        if abs(kpss_test_stats) >= abs(kpss_output['Critical Value 1%']):
            kpss_result = '99%'
        elif abs(kpss_test_stats) >= abs(kpss_output['Critical Value 2.5%']) and abs(kpss_test_stats) < abs(kpss_output[f'Critical Value 1%']):
            kpss_result = '97.5%'
        elif abs(kpss_test_stats) >= abs(kpss_output['Critical Value 5%']) and abs(kpss_test_stats) < abs(kpss_output[f'Critical Value 2.5%']):
            kpss_result = '95%'
        elif abs(kpss_test_stats) >= abs(kpss_output['Critical Value 10%']) and abs(kpss_test_stats) < abs(kpss_output[f'Critical Value 5%']):
            kpss_result = '90%'
    
        results_dct[col]['Differencing Stationary'] = adf_result
        results_dct[col]['Trending Stationary'] = kpss_result
        
    return results_dct


def detect_anomalies(squared_errors, n=1):
    threshold = np.mean(squared_errors) + n*np.std(squared_errors)
    detections = (squared_errors >= threshold).astype(int)
    
    return threshold, detections
mul_ts_df['Humidity'] = mul_ts_df['Humidity'].diff()
mul_ts_df['HumidityRatio'] = mul_ts_df['HumidityRatio'].diff()
mul_ts_df = mul_ts_df.dropna()
print(mul_ts_df.shape)

multi_ts_stationary = test_stationarity_multi_ts(mul_ts_df)
for k, v in multi_ts_stationary.items():
    print(k)
    print(v)
    print()<jupyter_output>(8142, 5)
Temperature
{'Differencing Stationary': '90%', 'Trending Stationary': '99%'}

Humidity
{'Differencing Stationary': '99%', 'Trending Stationary': '97.5%'}

Light
{'Differencing Stationary': '95%', 'Trending Stationary': '99%'}

CO2
{'Differencing Stationary': '99%', 'Trending Stationary': '99%'}

HumidityRatio
{'Differencing Stationary': '99%', 'Trending Stationary': '97.5%'}<jupyter_text>VAR to Detect Anomalies* The way it detects anomalies is to find observations with residuals above a threshold<jupyter_code># select better model order
max_lag = 20
var_model = VAR(mul_ts_df)
lag_results = var_model.select_order(max_lag)
selected_lag = lag_results.aic
print(f'Selected VAR order is {selected_lag}')

lag_results.summary()
model_fitted = var_model.fit(selected_lag)

# durbin_watson test to check whether there is any leftover pattern in the residuals, closer to 2, the better
dw_scores = durbin_watson(model_fitted.resid)

for col, dw in zip(mul_ts_df.columns, dw_scores):
    print(f'{col}: {dw}')
model_fitted.resid
squared_errors = model_fitted.resid.sum(axis=1)**2

threshold, detections = detect_anomalies(squared_errors, n=1)
detected_mul_ts_df = mul_ts_df.copy()
detected_mul_ts_df['anomaly_detection'] = detections
detected_mul_ts_df['Occupancy'] = occupancy
detected_mul_ts_df = detected_mul_ts_df.iloc[selected_lag:, :]
print(f'Threshold: {threshold}')

detected_mul_ts_df.head()
detected_mul_ts_df.loc[detected_mul_ts_df['anomaly_detection']==1].head()
# Check whether there's any anomaly pattern in different occupancy
no_occpupancy_df = detected_mul_ts_df.loc[detected_mul_ts_df['Occupancy']==0]
has_occpupancy_df = detected_mul_ts_df.loc[detected_mul_ts_df['Occupancy']==1]

print(no_occpupancy_df['anomaly_detection'].value_counts()/len(no_occpupancy_df))
print()
print(has_occpupancy_df['anomaly_detection'].value_counts()/len(has_occpupancy_df))<jupyter_output>0.0    0.995321
1.0    0.004679
Name: anomaly_detection, dtype: float64

0.0    0.991243
1.0    0.008757
Name: anomaly_detection, dtype: float64<jupyter_text>VECM to Detect Anomalies* About VECM: https://www.statsmodels.org/dev/generated/statsmodels.tsa.vector_ar.vecm.VECM.htmlstatsmodels.tsa.vector_ar.vecm.VECM<jupyter_code>k_ar_diff = 18
vecm_model = VECM(mul_ts_df, k_ar_diff=k_ar_diff)
vecm_model_fitted = vecm_model.fit()
vecm_dw_scores = durbin_watson(vecm_model_fitted.resid)

for col, dw in zip(mul_ts_df.columns, vecm_dw_scores):
    print(f'{col}: {dw}')
vecm_squared_errors = vecm_model_fitted.resid.sum(axis=1)**2

vecm_threshold, vecm_detections = detect_anomalies(vecm_squared_errors, n=1)
vecm_detected_mul_ts_df = mul_ts_df.iloc[k_ar_diff+1:, :]
vecm_detected_mul_ts_df['anomaly_detection'] = vecm_detections

print(f'Threshold: {threshold}')

vecm_detected_mul_ts_df.head()
compare_df = pd.merge(vecm_detected_mul_ts_df[['anomaly_detection']], detected_mul_ts_df[['anomaly_detection']], left_index=True, right_index=True)
print(len(compare_df))

compare_df.head()
compare_df.loc[(compare_df['anomaly_detection_x'] != compare_df['anomaly_detection_y'])]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13138440860215053, 'fraction_numerical': 0.021001344086021504, 'mean_word_length': 5.814539210074413, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 32, 'https://': 10, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '5325', 'license': 'MIT', 'n_tokens_mistral': 4674, 'n_tokens_neox': 4319, 'n_words': 744, 'path': 'code/yinyang/kats_experiments/kats_detect_outliers.ipynb', 'repo_name': 'lady-h-world/My_Garden', 'size': '11904'}