Pairs Trading with Kalman filter

Author

Lumpen95

1. 개요

공적분(cointegration)을 설명할 때 자주 인용되는 비유로 취한 사람과 그의 애완견 이야기가 있습니다. 구경꾼인 우리로서는 취한 사람과 애완견 모두 어디로 갈지 알 수 없습니다. 이러한 움직임은 랜덤 워크(random walk)의 전형적인 예입니다. 그렇지만 그 둘이 서로 줄로 연결되어 있다면 상황은 달라집니다. 애완견이 미친듯이 뛰어다닌다 해도 취한사람은 애완견과 일정한 거리를 유지하며 함께 이동할 것입니다. 두 시계열이 각각 비정상(non-stationary)임에도 불구하고 특정 선형 조합이 정상(stationary)을 이루는 관계, 즉 공적분 관계가 성립합니다.

공적분 관계에 있다는 것은 두 가격 시계열 간에 장기적 균형 관계가 유지된다는 의미이며, 이 균형이 일시적으로 붕괴될 때 이를 활용하여 수익을 추구하는 전략이 페어 트레이딩(pairs trading)입니다.

본 리포트에서는 KRX 헬스케어 지수에 포함된 개별 종목을 대상으로

  1. 공적분 관계가 가장 안정적으로 유지되는 페어를 통계적으로 산정하고,
  2. 칼만 필터(kalman filter)를 이용해 시간 가변(time-varying) \(\beta\) 를 추정한 후,
  3. 평균회귀(mean-reversion) 특성 기반의 트레이딩 전략을 구축 및 검증하는 것을 목표로 합니다.

2. 본론

2.1 종목 선택

  • 초기 투자 유니버스: KRX헬스케어 지수 내 개별 종목의 2023/01/02 기준과 2025/11/20 기준의 합집합

  • 전체 데이터 기간: 2023/01/02~2025/11/20 (3년)

  • 선정 방법

    1. 초기 투자 유니버스 종목으로 1:1 조합 데이터셋 생성
    2. 조합 별로 3년 데이터셋 Rolling k-fold 교차검증을 위한 split 생성
      • In-Sample(IS): 1 year
      • Out-Of-Sample(OSS): 3 months
    3. 다음 기준을 적용해서 조합 별로 점수를 부여하고 최고 점수를 받은 조합 선택
      1. Johansen 검정
        • IS / OSS 기간에 Johansen 검정을 실시하여 공적분 여부를 확인
        • 두개 기간 모두 공적분이 확인될 경우 1점, 아니면 0 점
      2. \(\beta\) Stability
        • Johansen 검정으로부터 IS \(\beta\) , OSS \(\beta\) 계산
        • \(\beta\) Stability = OSS \(\beta\) / IS \(\beta\)
        • 0.7 ~ 1.3 이내면 1점, 아니면 0점
      3. Half-life (평균 회귀 속도)
        • IS \(\beta\) 를 기준으로 spread 계산
          • \(spread_t = y_t - \beta x_t\)
        • spread에 대한 자기 회귀식 적용
          • \(\Delta spread_t = \alpha + \theta spread_{t-1} + \epsilon_t\)
        • \(\text{Half-life} = -\frac{ln(2)}{ln(1 + \theta)}\)
        • 50일 미만이면 1점, 100일 미만이면 0.5점, 이외 0점
    4. 총점은 0~3점이며, 최종적으로 최고 점수를 받은 1개 조합을 분석 대상으로 선정한다.

2.1.1 패키지 로드

Code 
# package load
library(quantmod)
library(tidymodels)
library(tidyverse)
library(timetk)
library(tseries)
library(urca)
library(tidyquant)
library(readxl)
library(readr)

options(scipen = 999)

2.1.2 KRX 헬스케어 데이터 로드

Code 
bio_2023 <- read_xlsx("KRX_bio_20230102.xlsx")
bio_2025 <- read_xlsx("KRX_bio_20251120.xlsx")

bio <- full_join(bio_2023, bio_2025, by = "종목코드") |> 
  mutate(name = ifelse(is.na(종목명.x), 종목명.y, 종목명.x)) |> 
  rename(code = 종목코드) |> 
  select(code, name)

bio
# A tibble: 102 × 2
   code   name            
   <chr>  <chr>           
 1 207940 삼성바이오로직스
 2 068270 셀트리온        
 3 091990 셀트리온헬스케어
 4 302440 SK바이오사이언스
 5 326030 SK바이오팜      
 6 000100 유한양행        
 7 128940 한미약품        
 8 137310 에스디바이오센서
 9 068760 셀트리온제약    
10 008930 한미사이언스    
# ℹ 92 more rows
  • 개별 종목 가격 테이블 생성
Code 
# 시작일, 종료일 설정
std <- as.Date("2023-01-02")
end <- as.Date("2025-11-20")

# 가격 호출 함수 정의
get_prc <- function(code) {
  
  ks_code <- paste0(code, ".KS")
  prc <- tryCatch(
    {
    getSymbols(ks_code, from = std, to = end, auto.assign = F, warnings = F, verbose = F) |> Ad() |> na.omit()
    },
    error = function(e) {
      kq_code <- paste0(code, ".KQ")
      tryCatch(
        {
          getSymbols(kq_code, from = std, to = end, auto.assign = F, warnings = F) |> Ad() |> na.omit()
        },
        error = function(e2) {
          return(NA)
        }
      )
    }) 
  return(prc)
}

len <- 703

bio <- bio |> 
  mutate(data = map(code, ~ get_prc(.x))) |> # 가격 호출 함수
  dplyr::filter(map_lgl(data, ~inherits(.x, "xts"))) |>  # 상장폐지 코드 제외
  dplyr::filter(map_lgl(data, ~length(.x) == len)) |>  # 시계열 길이 동일하게 설정
  mutate(
    log_data = map(data, ~log(.x)) # log price 변환
    )

bio
# A tibble: 92 × 4
   code   name             data            log_data       
   <chr>  <chr>            <list>          <list>         
 1 207940 삼성바이오로직스 <xts [703 × 1]> <xts [703 × 1]>
 2 068270 셀트리온         <xts [703 × 1]> <xts [703 × 1]>
 3 302440 SK바이오사이언스 <xts [703 × 1]> <xts [703 × 1]>
 4 326030 SK바이오팜       <xts [703 × 1]> <xts [703 × 1]>
 5 000100 유한양행         <xts [703 × 1]> <xts [703 × 1]>
 6 128940 한미약품         <xts [703 × 1]> <xts [703 × 1]>
 7 137310 에스디바이오센서 <xts [703 × 1]> <xts [703 × 1]>
 8 068760 셀트리온제약     <xts [703 × 1]> <xts [703 × 1]>
 9 008930 한미사이언스     <xts [703 × 1]> <xts [703 × 1]>
10 196170 알테오젠         <xts [703 × 1]> <xts [703 × 1]>
# ℹ 82 more rows

2.1.3 Pair set 생성

Code 
pairs <- t(combn(bio$name, 2)) |> 
  as_tibble()

pairs
# A tibble: 4,186 × 2
   V1               V2              
   <chr>            <chr>           
 1 삼성바이오로직스 셀트리온        
 2 삼성바이오로직스 SK바이오사이언스
 3 삼성바이오로직스 SK바이오팜      
 4 삼성바이오로직스 유한양행        
 5 삼성바이오로직스 한미약품        
 6 삼성바이오로직스 에스디바이오센서
 7 삼성바이오로직스 셀트리온제약    
 8 삼성바이오로직스 한미사이언스    
 9 삼성바이오로직스 알테오젠        
10 삼성바이오로직스 대웅제약        
# ℹ 4,176 more rows

2.1.4 Score function 정의

  • `get_pairs_prc : 정의된 Pair set 으로 가격 데이터셋 병합
Code 
get_pairs_prc <- function(V1, V2, data) {
  data |> 
    dplyr::filter(name %in% c(V1, V2)) |> 
    pull(log_data) |> 
    purrr::reduce(merge) |> 
    `colnames<-`(c(V1, V2))
}

  • get_cajo : Johansen 검정 함수

    • urca::ca.jo : Johansen 검정 함수로 다변량 시계열이 장기적으로 균형관계를 가지는지를 판별하기 위한 VAR 기반 기법

    • 다변량 시계열을 VAR 모델로 적합한 뒤 VECM으로 변환 이후 \(\Pi\) 행렬의 고유값을 이용해 공적분 여부 판단(고유값이 0보다 크면 공적분 관계 존개 가능)

    • Trace , Eigen statistic으로 공적분 개수 r 결정

      • x : 다변량 시계열 데이터 행렬

      • type : 공적분 개수 r 판정하는 방식(Trace , Eigen)

      • ecdet VECM의 deterministic term 결정

        • none : 상수 없음

        • const : 상수 포함(평균이 0이 아니거나 drift가 존재하는 경우)

        • trend : 추세 포함

      • K : VAR 차수(VECM에서는 K-1 lag 사용)

      • spec

        • longrun : 표준 Johansen VECM 방식

        • transitory : 단기 조정행렬이 다르게 계산되는 비표준 형태로 금융 분야에서는 거의 사용되지 않는 방식

  • get_cajo_score : Johansen 검정 결과 점수 집계 함수

Code 
get_cajo <- function(df) {
  tryCatch(
    ca.jo(df[2:3], type = "trace", ecdet = "none", K = 2, spec = "longrun"),
    error = function(e) NULL
  )
}

get_cajo_score <- function(cajo) {
  
  if (is.null(cajo)) return(0)
  h0_t_stat <- cajo@teststat[2]
  h0_5p <- cajo@cval[2, 2]
  
  if (is.na(h0_t_stat)) return(0)
  if (h0_t_stat > h0_5p) 1 else 0
}
  • `get_beta_score : \(\beta\) statibility 점수 집계 함수
Code 
get_beta_score <- function(cajo_is, cajo_oss) {
  
  if (is.null(cajo_is) | is.null(cajo_oss)) return(0)
  
  beta_is <- cajo_is@V[2, 1] / cajo_is@V[1, 1]
  beta_oss <- cajo_oss@V[2, 1] / cajo_oss@V[1, 1]
  
  ratio <- beta_oss / beta_is
  
  score <- if (ratio >= 0.7 && ratio <= 1.3) 1 else 0
}
  • get_half_life_score : Half-life 점수 집계 함수
Code 
get_half_life_score <- function(is_df, cajo_is) {
  if (is.null(cajo_is)) return(0)
  beta_is <- cajo_is@V[2, 1] / cajo_is@V[1, 1]
  spread_is <- is_df[, 3] - beta_is * is_df[, 2]
  spread_is <- na.omit(spread_is)
  m <- lm(diff(spread_is) ~ spread_is[-length(spread_is)])
  phi <- coef(m)[2]
  half_life <- -log(2) / log(1 + phi)
  
  if (is.na(half_life) || half_life <= 0) return(0)
  if (half_life < 50) {
    score <- 1
  } else if (half_life < 100) {
    score <- 0.5
  } else {
    score <- 0
  }
}
  • get_score : 최종 점수 집계 함수
Code 
get_score <- function(cv) {
  
  IS <- cv |> 
    pull(splits) |> 
    map(~analysis(.x))
  OSS <- cv |> 
    pull(splits) |> 
    map(~assessment(.x))
  
  ## score table
  cajo_tbl <- tibble(
    IS = IS,
    OSS = OSS,
    IS_cajo = map(IS, ~get_cajo(.x)),
    OSS_cajo = map(OSS, ~get_cajo(.x)),
    
    ## 1. Johanson test score
    IS_cajo_score = map_dbl(IS_cajo, ~get_cajo_score(.x)),
    OSS_cajo_score = map_dbl(OSS_cajo, ~get_cajo_score(.x)),
    cajo_score = as.integer(IS_cajo_score == 1 & OSS_cajo_score == 1),
    
    ## 2. beta stability
    beta_score = map2_dbl(IS_cajo, OSS_cajo, ~get_beta_score(.x, .y)),
    
    ## 3. half life
    half_life_score = map2_dbl(IS, IS_cajo, ~get_half_life_score(.x, .y))
    
  )
  
  final_score <- cajo_tbl |> 
    select(cajo_score, beta_score, half_life_score) |> 
    colMeans() |> 
    sum()
}

2.1.5 최종 pairs 선택

정의한 함수들을 사용하여 최종 pairs를 선택하고 assessment 가 3개월인 k-fold split을 생성합니다.

Code 
pairs <- pairs |> 
  mutate(pairs_prc = map2(V1, V2, ~get_pairs_prc(.x, .y, bio)),
         pairs_prc = map(pairs_prc, ~fortify.zoo(.x)),
         cv = map(pairs_prc, ~timetk::time_series_cv(.x, 
                                                     date_var = Index,
                                                     initial = "1 years",
                                                     assess = "3 months",
                                                     skip = "3 months")))
pairs <- pairs |> 
  mutate(coint_score = map_dbl(cv, ~get_score(.x)))

best_pairs <- 
  pairs |> 
  slice_max(coint_score)

best_pairs
  • 한미약품과 한미사이언스가 1.71 점으로 최종 선택되었습니다.

먼저 고정 \(\beta\) 를 사용해 spread를 계산해보겠습니다. ca.jo 함수를 통해 계산이 가능합니다.

Code 
best_df <- best_pairs |> 
  pluck("pairs_prc", 1) |> 
  as_tibble() |> 
  mutate(across(contains("한미"), exp))

cajo <- ca.jo(best_df[, 2:3], type = "trace", ecdet = "none", K = 2, spec = "longrun")
fixed_beta <- cajo@V[2, 1]

fixed_beta_df <- best_df |> 
  mutate(
    fixed_spread = 한미약품 - fixed_beta * 한미사이언스
  )

# 시각화
fixed_beta_df |> 
  pivot_longer(-Index) |> 
  mutate(name = factor(name, levels = c("한미약품", "한미사이언스", "fixed_spread"))) |> 
  plot_time_series(.date_var = Index, 
                   .interactive = F, 
                   .value = value, 
                   .smooth = F, 
                   .facet_vars = name, 
                   .title = expression("Fixed " * beta), 
                   .color_var = name, 
                   .legend_show = F)

Code 
# 통계적 검정
fixed_beta_df |> 
  pull(fixed_spread) |> 
  adf.test()

    Augmented Dickey-Fuller Test

data:  pull(fixed_beta_df, fixed_spread)
Dickey-Fuller = -2.3276, Lag order = 8, p-value = 0.4396
alternative hypothesis: stationary

ca.jo 함수를 통해 계산된 \(\text{Fixed}~\beta\) 로 spread를 계산한 결과 다음과 같은 결론을 내릴 수 있다.

  • 검정 결과 랜덤워크는 아니다.

  • 페어 트레이딩에 적합한 mean-reverting 의 모습을 보인다.

  • 초중반 기간은 대체로 양호한 공적분 관계를 보이고 있었지만, 후반부에는 신호가 약해지고 있다.

  • \(\text{Fixed}~\beta\) 임에도 나름 괜찮은 공적분 관계를 보이고 있다.

  • \(\text{Time-varying}~\beta\) 를 사용할 경우 더 훌륭한 공적분 관계를 보일 것을 기대할 수 있다.

2.2 \(\text{Time-varying}~\beta\) estimation using kalman filter

2.2.1 kalman filter

칼만 필터는 대표적인 시계열 상태공간 모델(State-Space model)입니다. 상태공간 모델에서는 참 상태(true state)를 직접 측정할 수 없으며 오직 측정(observation)된 것으로부터 추론하는 것만 가능합니다. 즉, 상태(state)는 시간에 따라 변하며 그 상태로부터 생성된 관측값으로 상태를 추론할 수 있는 것입니다.

관측에 기반한 추정 작업은 다음과 같이 구분될 수 있습니다.

  • 필터링: 시간 \(t\) 의 상태에 대한 추정 갱신에 시간 \(t\) 의 측정 사용 - 현재

  • 예측: 시간 \(t\) 의 예상되는 상태에 대한 예측 생성에 시간 \(t-1\) 의 측정을 사용 - 과거

  • 평활화: 시간 \(t\) 의 참 상태 추정에 시간 \(t\) 및 \(t\) 의 전후 측정 사용 - 현재 + 과거

칼만 필터는 Noise가 있는 관측값으로부터 숨겨진 상태 변수를 실시간으로 추정하는 알고리즘으로 계산이 간단하고, 미래 예측이나 현재를 추정하는 데 과거 데이터의 저장이 필요 없다는 것이 장점입니다.

시장 환경이 변화함에 따라 자산 간 민감도에 해당하는 \(\beta\) 역시 시간에 변합니다. 따라서 백테스팅 뿐만 아니라 실제 전략 운용 단계에서도 고정된 \(\beta\) 를 사용하면 시장 변동을 충분히 반영하지 못하는 문제가 발생합니다. 이에 계산 효율적이면서도 노이즈에 강한 칼만 필터를 활용하여 시변 \(\beta\) 를 추정하고 이를 전략에 반영하고자 합니다.

칼만 필터의 계산 과정은 다음과 같습니다. 칼만 필터는 많은 양의 추적이 필요하기에 순환적인 반복성을 띕니다.

  1. 모델 정의
    • 칼만 필터는 사전에 다음 2가지 모델의 정의가 필요합니다.
    • State transition model
      • 상태 변수가 다음 상태 변수로 변하는 과정을 설명하는 모델
      • \(X_t = FX_{t-1} + w_t,~w_t \sim N(0, Q)\)
    • Observation model
      • \(Y_t = HX_t + v_t, v_t \sim N(0, R)\)
  2. 상태 예측
    • \(\hat{X_{t|t-1}} = F \hat{X_{t-1|t-1}}\)
  3. 상태 공분산 예측
    • \(P_{t|t-1} = FP_{t-1|t-1}F^T + Q\)
  4. 관측 예측 오차
    • \(v_t = Y_t - H\hat{X_{t|t-1}}\)
  5. 관측 예측 오차 공분산
    • \(S_t = HP_{t|t-1}H^T + R\)
  6. 칼만 이득(kalman gain)
    • \(K_t = \frac{P_{t|t-1}H^T}{S_t}\)
  7. 상태 갱신(Update)
    • \(\hat{X_{t|t}} = \hat{X_{t|t-1}} + K_tv_t\)
  8. 상태 공분산 갱신
    • \(P_{t|t} = (I - K_tH)P_{t|t-1}\)
  9. 다시 1번으로 돌아가서 추정과 갱신 반복

2.2.2 \(\text{Time-varying}~\beta\) estimation

칼만필터를 사용해 \(\beta\) 를 추정하는 경우 일반적으로 다음과 같이 랜덤워크를 가정하여 상태모델과 관측모델을 정의합니다.

  • \(\beta_t = \beta_{t-1} + w_t,~w_t \text{~}N(0, Q)\)

  • \(y_t = x_t\beta_t + v_t, v_t \text{~} N(0, R)\)

칼만필터와 관련된 다양한 R 패키지가 존재하지만 수동으로 구현해보려고 합니다.

  • \(y\) : 한미약품(종속변수)

  • \(x\) : 한미사이언스(독립변수)

Code 
## 시계열 변수 정의
y <- best_df$한미약품
x <- best_df$한미사이언스
n <- nrow(best_df)

## kalman filter 함수 정의
# 상태
beta <- rep(0, n)
# 상태 공분산
P <- rep(0, n)
# 관측 예측 오차
pred_error <- rep(0, n)
# 관측 예측 오차 공분산
Q <- rep(0, n)

# # observation noise
sigma_epsilon <- 1e-2 
ratio <- 0.05
# state noise
sigma_eta <- sqrt(ratio) * sigma_epsilon 
# 공분산 최소값 제한
P_min <- 1e-4

beta[1] <- 0
P[1] <- 1
Q[1] <- 0

for (t in 2:n) {
  # Prediction
  beta_pred <- beta[t-1]
  P_pred <- P[t-1] + sigma_eta^2
  pred_error[t] <- y[t] - x[t] * beta_pred
  Q[t] <- x[t]^2 * P_pred + sigma_epsilon^2
  
  # Kalman gain
  K <- (P_pred * x[t]) / Q[t]
  
  # Update
  beta[t] <- beta_pred + K * pred_error[t]
  P[t] <- (1 - K * x[t]) * P_pred
  
  if (P[t] < P_min) P[t] <- P_min
}
Code 
# 시변 베타 시계열 테이블 생성
kalman_df <- best_df |> 
  mutate(
    tv_beta = beta, # time-varying beta
    tv_spread = 한미약품 - tv_beta * 한미사이언스,
    tv_error = pred_error,
    tv_sqrtQ = sqrt(Q)
  ) |> 
  tail(-2) # 초기 beta는 0이기에 제외

# 시각화
kalman_df |> 
  select(-tv_beta, -tv_error, - tv_sqrtQ) |> 
  pivot_longer(-Index) |> 
  mutate(name = factor(name, levels = c("한미약품", "한미사이언스", "tv_spread"))) |> 
  plot_time_series(.date_var = Index, 
                   .interactive = F, 
                   .value = value, 
                   .smooth = F, 
                   .facet_vars = name, 
                   .title = expression("Time-Varying " * beta * " using Kalman Filter"), 
                   .color_var = name, 
                   .legend_show = F)

Code 
# 통계 검정
kalman_df |> 
  pull(tv_spread) |> 
  adf.test()

    Augmented Dickey-Fuller Test

data:  pull(kalman_df, tv_spread)
Dickey-Fuller = -9.6925, Lag order = 8, p-value = 0.01
alternative hypothesis: stationary

  • tv_spread는 fixed_spread 와 달리 랜덤워크 모습을 보이며 잔차 그래프와 비슷해 보인다.

  • 페어 트레이딩 시 매매 시그널 생성을 위해서는 스프레드의 정상성이 필수적이다. 그러므로 랜덤워크인 tv_spread가 페어 트레이딩에 더 적합한 것을 알 수 있다.

2.3 Trading Signal

추정한 \(\beta\) 를 사용해 매매 신호 규칙을 정한다. 일반적으로는 계산된 tv_spread의 rolling z-score를 통해 신호를 생성하지만 이번에는 안정성이 더 높다고 평가되는 \(\beta\) 를 추정할 때 계산한 관측 예측 오차(\(v_t\))를 사용해 신호를 생성할 것이다.

규칙은 다음과 같다.

\[ \begin{split} z_t &= \frac{v_t - \mu_t}{\sigma_t} \\ \mu_t &= \frac{1}{L}\sum^t_{i=t-L+1} v_i \\ \sigma_t &= \sqrt{\frac{1}{L-1}\sum^t_{i=t-L+1} (v_i - \mu_t)^2} \\ \text{signal}_t &= \begin{cases} -1 , & z_t > z_{\text{entry}} \\ +1 , & z_t < -z_{\text{entry}} \end{cases} \end{split} \]

  • 예측 오차(\(v_t\)) 를 대상으로 60일 rolling-window 방식으로 z-score(\(z_t\)) 계산

  • 매매 시그널 생성에 사용하기 위한 optimal threshold(\(z_{entry}\)) 는 Palomar 교수님의 Non-parametric 방식으로 추정

    • Pairs Trading - Daniel P. Palomar

  • Non-parametric optimal threshold 추정 과정은 다음과 같습니다.

    1. Observed sample path

      • \(z_1, z_2, ..., z_T\)

    2. Discretized threshold values

      • \(s_o \in \{s_{o1}, s_{o2}, ..., s_{oJ}\}\)

    3. Emperical trading frequency :

      • \[ \bar{f_j} = \frac{\sum_{t=1}^T1\{Z_t>s{oj\}}}{T} \]

    4. Smooth the trading frequency function by regularization

      • \[ \text{minimize}_f \quad \sum_{j=1}^J (\bar{f}_j - f_j)^2 + \lambda \sum_{j=1}^{J-1} (f_j - f_{j+1})^2 \]

      • 위 식은 다음과 같이 바꿔쓸 수 있습니다.

      • \[ \text{minimize}_f \quad \lVert\bar{f}-f\rVert^2_2 + \lambda\lVert Df \rVert ^2_2 \\ \]

      • \[ D = \begin{bmatrix} 1 & -1 & & \\ & 1 & -1 & \\ & & \ddots & \ddots & \\ & & & 1 & -1 \end{bmatrix} \in \mathbb{R}^{(J-1) \times J} \]

      • \[ f^{\star} = (I + \lambda D^T D)^{-1\bar{f}} \]

    5. The optimal threshold is the one maximizes the total profit

      • \[ S_o^{\star} = \underset{s_{oj} \in \{s_{o1}, s_{o2}, \ldots, s_{oJ}\}}{\arg\max} \{s_{oj}f_j\} \]
Code 
# 예측 오차
err_vec <- kalman_df$tv_error

# 최적 임계값 생성 함수
make_s0 <- function(err_vec, window, lambda = 10, grid_n = 200) {
  
  m <- length(err_vec)
  z_t <- rep(NA, m)
  s_0 <- rep(NA, m)
  
  for (i in 1:m) {
    
    if (i >= window) {
      
      # window 기간 정의
      window_ts <- err_vec[(i - window + 1):i]
      
      # z-score
      window_z <- (window_ts - mean(window_ts)) / sd(window_ts)
      z_t[i] <- window_z[window]
      
      # optimal threshold
      min_s <- quantile(abs(window_z), probs = 0.75)
      max_s <- quantile(abs(window_z), probs = 0.95)
      
      # grid
      s0_grid <- seq(min_s, max_s, length.out = grid_n)
      J <- length(s0_grid)
      
      # empirical trading frequency
      f_emp <- map_dbl(s0_grid, ~ mean(abs(window_z) > .x))
      D <- diff(diag(J), differences = 1)
      I <- diag(J)
      
      # smoothing
      f_smooth <- solve(I + lambda * t(D) %*% D, f_emp)
      
      # profit
      profit <- s0_grid * f_smooth
      
      optimal_idx <- which.max(profit)
      
      s_0[i] <- s0_grid[optimal_idx]
    }
  }
      return(
        tibble(
          z_t = z_t,
          s_0 = s_0,
          neg_s_0 = -s_0
        )
      )
}

signal_tbl <- make_s0(err_vec, window = 60)
signal_xts <- xts(signal_tbl, order.by = kalman_df$Index)

# signal 분류
signal_xts$signal <- ifelse(
    (signal_xts$z_t > signal_xts$s_0 & lag.xts(signal_xts$z_t, 1) < lag.xts(signal_xts$s_0, 1)), 
    -1,
      ifelse(
        (signal_xts$z_t < signal_xts$neg_s_0 & lag.xts(signal_xts$z_t, 1) > lag.xts(signal_xts$neg_s_0, 1)),
              1, 0)
    )

sig <- signal_xts$signal
sig[sig == 0] <- NA
sig <- na.locf(sig)
sig <- diff(sig) / 2

sig <- sig |> 
  as.numeric() |> 
  na_if(0)

2.4 Backtesting

현재까지 생성된 시그널은 다음과 같은 규칙을 가지고 있습니다.

  • z-score가 최적 임계값보다 작았지만 커진 경우 => -1 (short spread - short y, long x)
  • z-score가 최적 임계값보다 컸지만 작아진 경우 => +1 (long spread - long y, short x)
  • 나머지는 모두 NA로 채워진 상태

이제부터 우리는 다음과 같은 사항을 고려해 시그널 생성을 마무리 해야합니다.

  • Exit 시점 추가(sig = 0)

    • 포지션 매수와 매도만 고려할 경우 주식 간 공적분 관계가 무너질 경우 손실이 무한대로 커질 수 있습니다.

    • 이를 방지하기 위해 예측오차(tv_error)를 활용한 안전장치를 추가할 것입니다.

    • \(tv\_error > tv\_error\_rolling\_median + tv\_error\_rolling\_sd \times 2\)

    • 위 조건을 만족하는 경우 포지션을 청산합니다.

  • 시그널이 전환되지 않는 한 기존 시그널을 유지하도록 설정합니다.

  • 사전관찰을 피하기 위해 t 일에 생성된 시그널을 t + 1 일에 사용합니다.

위 고려사항을 모두 포함한 시그널 및 누적 손익은 다음과 같습니다.

Code 
PnL_df <- kalman_df |> 
  mutate(sig = lag(sig)) |> # 사전 관찰 방지
  tq_mutate(
    tv_error, mutate_fun = runSD, col_rename = "tv_error_rolling_sd", n = 60
  ) |> 
  tq_mutate(
    tv_error, mutate_fun = runMedian, col_rename = "tv_error_rolling_median", n = 60
  ) |> 
  mutate(
    sig = ifelse(
      (abs(tv_error) > tv_error_rolling_median + tv_error_rolling_sd * 2), 0, sig
    ),
    sig = na.locf(sig, na.rm = F) |> 
      replace_na(0),
    dt_y = c(0, diff(한미약품)),
    dt_x = c(0, diff(한미사이언스)),
    PnL = sig * (dt_y - tv_beta * dt_x),
    cumPnL = cumsum(PnL)
  )

PnL_df <- PnL_df |> 
  select(Index, 한미약품, 한미사이언스, tv_beta, tv_error, sig, cumPnL)
Code 
PnL_df |> 
  select(-sig, -cumPnL) |> 
  pivot_longer(-Index) |>
  mutate(
    name = factor(name, 
                  levels = c("한미약품", 
                             "한미사이언스", 
                             "tv_beta", 
                             "tv_error"), 
                  labels = c("한미약품 (Y)", 
                             "한미사이언스 (X)", 
                             "Time-Varying Beta", 
                             "kalman error"))
    ) |> 
  ggplot(aes(Index, value, col = name)) +
  geom_line(show.legend = F) +
  facet_wrap(~name, ncol = 1, scales = "free") +
  scale_y_continuous(labels = scales::comma) +
  theme_tq()

Code 
PnL_df %>%
  ggplot(aes(Index)) +
  
  geom_line(aes(y = cumPnL, color = "cumPnL"), linewidth = .6) +
  
  geom_segment(
    data = subset(PnL_df, sig == 0),
    aes(x = Index - 1, xend = Index, y = cumPnL, yend = cumPnL, color = "Neutral (sig = 0)"),
    , linewidth = 1.3, lineend = "round") +
  
  geom_line(aes(y = tv_error, color = "kalman error"), alpha = .4, linewidth = .5) +
  
  geom_hline(yintercept = 0, color = "grey50", linewidth = .4) +
  
  scale_y_continuous(
    name = "cumPnL",
    labels = scales::comma
  ) +
  
  scale_color_manual(
    values = c(
      "cumPnL" = "black",
      "Neutral (sig = 0)" = "red",
      "kalman error" = "steelblue"
    )
  ) +
  
  labs(
    title = "Pairs Trading PnL & Time-Varying Beta error",
    x = NULL, 
    color = NULL
  ) +
  
  geom_rect(
      aes(xmin = as.Date("2023-10-01"), xmax = as.Date("2024-05-01"), ymin = -100000, ymax = 70000), 
      fill = NA, 
      color = "black", 
      linewidth = 1
  ) +
  
  annotate(
    "label",
    x = as.Date("2023-10-20"), y = 77000,
    label = "1",
    size = 6,
    fontface = "bold"
  ) +
  
  geom_rect(
    aes(xmin = as.Date("2024-08-01"), xmax = as.Date("2025-01-01"), ymin = -50000, ymax = 120000), 
    fill = NA, 
    color = "black", 
    linewidth = 1
  ) +
  
  annotate(
    "label",
    x = as.Date("2024-08-20"), y = 127000,
    label = "2",
    size = 6,
    fontface = "bold"
  ) +
  
  theme_tq() +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    legend.position = "inside",
    legend.justification = c(0.03, 0.98),
    legend.background = element_rect(color = "black"),
    legend.text = element_text(size = 12)
  )

  • 2023년 초기 약 6개월은 전략 추정(rolling kalman filter 기반 \(\beta\) 안정화) 를 위한 관찰 기간으로 의도적으로 포지션을 취하지 않았습니다.

  • 1번 변동성 구간 이전까지는 스프레드(kalman error)가 장기간 평균 부근에서 움직이면서 전략 특성 상 수익 기회가 제한적이었습니다. 그러나 1번 구간에서 스프레드가 단기간 크게 이탈하면서 잔차 분산이 확대되었고, 이후 정상구간으로의 평균회귀가 강하게 나타나 높은 수익이 발생했습니다. 이는 공적분 기반 페어 트레이딩의 전형적 payoff 구조가 가장 잘 작동한 구간이라고 볼 수 있습니다.

  • 2번 변동성 구간 역시 1번 구간과 마찬가지로 스프레드의 급격한 변동이 관측되었고 일시적 붕괴로 판단되어 포지션이 청산되었습니다. 결과적으로 변동성 스파이크에 의한 잠재적 손실을 효과적으로 차단했습니다.

  • 1번 구간과 달리 2번 구간에서 매매 제한이 발생한 이유는 변동성의 절대적인 크기보다 \(\beta\) 추정 안정성이 크게 훼손되었기 때문으로 판단됩니다.

  • 2025년 이후로는 스프레드가 비교적 안정적인 범위에서 움직이고 있습니다.

3. 결론 : 상태

예측하지 마라 - 벤자민 그레이엄

수많은 투자 대가들이 입을 모아 시장을 예측하지 말라고 말합니다. 시장 자체를 예측하려는 시도는 의미가 없으며, 우리가 할 수 있는 것은 ‘현재 상태를 정확히 읽고 이에 대응하는 것’뿐이라는 점입니다.

칼만 필터를 적용한 페어 트레이딩 전략은 이러한 철학과 가장 잘 맞는 기법입니다. 칼만 필터는 시시각각 변하는 스프레드의 공적분 관계와 베타의 추정 오차를 실시간으로 업데이트하며, 시장에 내재된 상태(state)를 추정하고 그 상태에 맞춰 전략을 수정할 수 있게 합니다.

이 전략은 미래를 예측하려 하지 않습니다. 대신, 공적분 관계가 유지되는지, 잔차가 정상성을 벗어났는지, 베타 추정이 안정적인지와 같은 현재 시장의 구조적 정보를 반영해 포지션을 확대하거나 축소하며 위험을 통제합니다.

우리가 확실히 알고 있는 한 가지는 미래는 예측할 수 없다는 것이며, 우리가 할 수 있는 단 하나는 현재의 상태를 정확히 읽고 변화가 감지되면 즉시 대응하는 것입니다.

칼만 필터 기반 페어 트레이딩은 이러한 상태 기반 대응 철학을 구현한 전략이며, 복잡한 시장 환경 속에서도 일관된 규칙과 적응성을 제공한다는 점에서 의미가 있습니다.