# ==========================================
# 📊 AirPassengers 시계열 데이터 전처리
# ==========================================

import pandas as pd
import matplotlib.pyplot as plt

# 1️⃣ CSV 파일 불러오기
# Kaggle에서 받은 AirPassengers.csv 파일을 불러온다.
# 데이터에는 'Month'와 'Passengers' 두 개의 열이 포함되어 있다.
df = pd.read_csv("AirPassengers.csv")

# 2️⃣ 데이터 상위 5개 확인
print("원본 데이터 (상위 5행):")
print(df.head())

# 3️⃣ 'Month' 열을 날짜(datetime) 형식으로 변환
# 문자열 형식의 날짜를 datetime 형식으로 바꿔야
# 시계열 연산(예: 차분, 이동평균 등)을 수행 가능
df["Month"] = pd.to_datetime(df["Month"])

# 4️⃣ 'Month'를 인덱스로 설정하고, 월 단위(Month Start) 시계열로 지정
# asfreq("MS")는 "Month Start" 빈도로 시계열 주기 설정
df = df.set_index("Month").asfreq("MS")

# 5️⃣ 변환된 데이터 확인
print("\n시계열 인덱스 적용 후 데이터:")
print(df.head())

# 6️⃣ 데이터 기본 통계 요약
print("\n📈 데이터 통계 요약:")
print(df.describe())

# 7️⃣ 결측치 확인
print("\n결측치 확인:")
print(df.isnull().sum())

# 8️⃣ 시각화: 월별 여객 수 변화 추이
plt.figure(figsize=(10, 5))
plt.plot(df.index, df["#Passengers"], color="royalblue", linewidth=2)
plt.title("Monthly Number of Air Passengers (1949–1960)", fontsize=13)
plt.xlabel("Year")
plt.ylabel("Passengers")
plt.grid(True, linestyle="--", alpha=0.5)
plt.tight_layout()
plt.show()

# 9️⃣ 간단한 통계 해석 출력
print("\n📊 해석:")
print("• 데이터 기간:", df.index.min().strftime("%Y-%m"), " ~ ", df.index.max().strftime("%Y-%m"))
print("• 전체 관측치 수:", len(df))
print("• 월별 평균 여객 수:", round(df['#Passengers'].mean(), 2))
print("• 최대 여객 수:", df['#Passengers'].max())
print("• 최소 여객 수:", df['#Passengers'].min())

원본 데이터 (상위 5행):
     Month  #Passengers
0  1949-01          112
1  1949-02          118
2  1949-03          132
3  1949-04          129
4  1949-05          121
---------------------------

시계열 인덱스 적용 후 데이터:
            #Passengers
Month
1949-01-01          112
1949-02-01          118
1949-03-01          132
1949-04-01          129
1949-05-01          121
---------------------------

📈 데이터 통계 요약:
       #Passengers
count   144.000000
mean    280.298611
std     119.966317
min     104.000000
25%     180.000000
50%     265.500000
75%     360.500000
max     622.000000
---------------------------

결측치 확인:
#Passengers    0
dtype: int64
---------------------------

📊 해석:
• 데이터 기간: 1949-01  ~  1960-12
• 전체 관측치 수: 144
• 월별 평균 여객 수: 280.3
• 최대 여객 수: 622
• 최소 여객 수: 104
---------------------------

pip install openpyxl

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

# 1) 엑셀 파일 로드 (openpyxl 패키지가 필요함)
# - AirQualityUCI.xlsx 파일을 불러옴
# - sheet_name=0 : 첫 번째 시트를 읽음
# - engine="openpyxl" : .xlsx 파일을 처리하는 데 필요한 엔진 지정
df = pd.read_excel(
    "AirQualityUCI.xlsx",
    sheet_name=0,
    engine="openpyxl"
)

# 2) Date/Time 정규화
#   - Date: 어떤 타입이든 datetime으로 변환하여 날짜 형식을 통일
#   - dayfirst=True : 날짜가 'DD/MM/YYYY' 형식일 때 올바르게 인식되도록 지정
#   - errors="coerce" : 변환 불가능한 값은 NaT(Not a Time)으로 처리하여 오류 방지
date_ser = pd.to_datetime(df["Date"], dayfirst=True, errors="coerce")

# datetime 타입의 날짜를 문자열 'YYYY-MM-DD' 형식으로 변환
# 문자열로 만들어 두면 나중에 Time과 안전하게 결합할 수 있음
date_str = date_ser.dt.strftime("%Y-%m-%d")

# Time 컬럼 정규화
#   - Time 값이 HH.MM.SS 형식으로 되어 있으므로, 점(.)을 콜론(:)으로 바꿈
#   - astype(str) : 숫자나 결측값 등 모든 값을 문자열로 통일
#   - str.replace(".", ":", n=2, regex=False) : 앞에서 두 번만 '.'을 ':'로 교체 (HH.MM.SS → HH:MM:SS)
#   - str.strip() : 앞뒤 공백을 제거
time_str = (
    df["Time"]
    .astype(str)
    .str.replace(".", ":", n=2, regex=False)
    .str.strip()
)

# Time 문자열 중 길이가 0인 항목(즉, 비어 있는 값)을 찾아서
# "00:00:00"으로 채워주는 코드임.
#   -> Time 정보가 비어 있으면 자정(00시)으로 대체하여 Datetime 생성 시 오류를 방지
time_str = time_str.where(time_str.str.len() > 0, "00:00:00")

# 3) 최종 Datetime 생성 (문자열 결합 후 파싱)
df["Datetime"] = pd.to_datetime(date_str + " " + time_str, errors="coerce")

# 4) 인덱스/주기 설정
df = df.set_index("Datetime").sort_index().asfreq("H")  # 시간 단위(Hourly)

# 5) 센서 오류값 -200을 결측치로 처리 (UCI Air Quality 관례)
num_cols = df.select_dtypes(include="number").columns
df[num_cols] = df[num_cols].replace(-200, np.nan)

# 6) 확인
print(df.tail())
print("\n결측치 개수:")
print(df.isnull().sum())

# 7) 간단 시각화: CO(GT) 농도 시계열 그래프
if "CO(GT)" in df.columns:
    ax = df["CO(GT)"].plot(figsize=(12, 4), title="CO(GT) Concentration Over Time")
    ax.set_xlabel("Datetime")
    ax.set_ylabel("CO (mg/m^3)")
    plt.tight_layout()
    plt.show()
else:
    print("Column 'CO(GT)' not found in the DataFrame.")

[5 rows x 15 columns]

결측치 개수:
Date                0
Time                0
CO(GT)           1683
PT08.S1(CO)       366
NMHC(GT)         8443
C6H6(GT)          366
PT08.S2(NMHC)     366
NOx(GT)          1639
PT08.S3(NOx)      366
NO2(GT)          1642
PT08.S4(NO2)      366
PT08.S5(O3)       366
T                 366
RH                366
AH                366
dtype: int64

Date                0
Time                0
CO(GT)           1683
PT08.S1(CO)       366
NMHC(GT)         8443
C6H6(GT)          366
PT08.S2(NMHC)     366
NOx(GT)          1639
PT08.S3(NOx)      366
NO2(GT)          1642
PT08.S4(NO2)      366
PT08.S5(O3)       366
T                 366
RH                366
AH                366

# 결측치 보간 예시
df_interp = df.copy()

# (1) 선형 보간법으로 연속 구간 채우기
df_interp = df_interp.interpolate(method="linear", limit_direction="both")

# (2) 여전히 남은 결측치는 ffill, bfill로 보완
df_interp = df_interp.fillna(method="ffill").fillna(method="bfill")

# 보정 후 결측치 재확인
print("보정 후 결측치 합계:", df_interp.isnull().sum().sum())

보정 후 결측치 합계: 0

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

# 1️⃣ 엑셀 파일 로드
df = pd.read_excel("AirQualityUCI.xlsx", sheet_name=0, engine="openpyxl")

# 2️⃣ 시간 정보 정규화
df["Time"] = (
    df["Time"].astype(str)
    .str.replace(".", ":", n=2, regex=False)
    .str.strip()
    .where(lambda s: s.str.len() > 0, "00:00:00")
)
date_ser = pd.to_datetime(df["Date"], dayfirst=True, errors="coerce")
date_str = date_ser.dt.strftime("%Y-%m-%d")
df["Datetime"] = pd.to_datetime(date_str + " " + df["Time"], errors="coerce")

# 3️⃣ Datetime 인덱스 설정 및 시간 주기 지정
df = df.set_index("Datetime").sort_index().asfreq("h")

# 4️⃣ 센서 오류값(-200)을 NaN으로 변환
num_cols = df.select_dtypes(include=np.number).columns
df[num_cols] = df[num_cols].replace(-200, np.nan)
df[num_cols] = df[num_cols].infer_objects(copy=False)

# 5️⃣ 결측치 보간 (선형보간 + ffill + bfill)
df_interp = df[num_cols].interpolate(
    method="linear",
    limit_direction="both" # "forward" + "backward" 둘 다 고려해서 양쪽 값의 평균 경로로 보간하는 방식
).ffill().bfill()

# 6️⃣ 시각화 대상 변수 선택
target_col = "CO(GT)"
before = df[target_col]
after = df_interp[target_col]

# 7️⃣ 결측이 있었던 구간 찾기 (시각적 강조용)
mask_nan = before.isna()
if mask_nan.sum() > 0:
    first_nan = mask_nan[mask_nan].index.min()
    last_nan = mask_nan[mask_nan].index.max()
    start_zoom = first_nan - pd.Timedelta(days=2)
    end_zoom = last_nan + pd.Timedelta(days=2)
else:
    start_zoom = df.index.min()
    end_zoom = start_zoom + pd.Timedelta(days=7)

# 8️⃣ 상하 2단 subplot 시각화
fig, axes = plt.subplots(2, 1, figsize=(14, 7), sharex=True)

# --- 상단: 보정 전 ---
axes[0].plot(before.index, before, color="gray", alpha=0.6, label="Before Correction")
axes[0].set_title(f"{target_col} - Before Correction")
axes[0].set_ylabel("CO (mg/m³)")
axes[0].grid(alpha=0.3)
axes[0].legend()

# --- 하단: 보정 후 ---
axes[1].plot(after.index, after, color="blue", linewidth=1.5, label="After Interpolation + Fill")
axes[1].set_title(f"{target_col} - After Interpolation + Fill")
axes[1].set_xlabel("Datetime")
axes[1].set_ylabel("CO (mg/m³)")
axes[1].grid(alpha=0.3)
axes[1].legend()

# 확대 영역 설정
axes[0].set_xlim(start_zoom, end_zoom)
axes[1].set_xlim(start_zoom, end_zoom)

plt.tight_layout()
plt.show()

                          Date      Time  CO(GT)  PT08.S1(CO)  ...  PT08.S5(O3)       T      RH
  AH
Datetime                                                       ...

2005-04-04 10:00:00 2005-04-04  10:00:00     3.1      1314.25  ...      1728.50  21.850  29.250  0.756824
1864
2005-04-04 12:00:00 2005-04-04  12:00:00     2.4      1142.00  ...      1092.00  26.900  18.350  0.640649
2005-04-04 13:00:00 2005-04-04  13:00:00     2.1      1002.50  ...       769.75  28.325  13.550  0.513866
2005-04-04 14:00:00 2005-04-04  14:00:00     2.2      1070.75  ...       816.00  28.500  13.125  0.502804

[5 rows x 15 columns]

보정 후 전체 결측치 수: 0

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

# 1. 데이터 로드 및 준비
df = pd.read_excel("AirQualityUCI.xlsx", sheet_name=0, engine="openpyxl")
df["Time"] = df["Time"].astype(str).str.replace(".", ":", n=2, regex=False).str.strip()
date_ser = pd.to_datetime(df["Date"], dayfirst=True, errors="coerce")
date_str = date_ser.dt.strftime("%Y-%m-%d")
df["Datetime"] = pd.to_datetime(date_str + " " + df["Time"], errors="coerce")
df = df.set_index("Datetime").sort_index().asfreq("h")

# 2. 숫자형 컬럼만 추출 및 센서 오류값(-200) 제거
num_cols = df.select_dtypes(include=np.number).columns
df[num_cols] = df[num_cols].replace(-200, np.nan)

# 3. 대상 컬럼 선택
target_col = "CO(GT)"
series = df[target_col].dropna()

# 4. IQR 계산
Q1 = series.quantile(0.25)
Q3 = series.quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR

# 5. 이상치 탐지
outliers = series[(series < lower_bound) | (series > upper_bound)]

print(f"이상치 개수: {len(outliers)}")
print(f"상한: {upper_bound:.2f}, 하한: {lower_bound:.2f}")

# 6. 시각화
plt.figure(figsize=(12,5))
plt.plot(series.index, series, label="Original", color="gray")
plt.scatter(outliers.index, outliers, color="red", label="Outliers")
plt.title("Outlier Detection using IQR Method")
plt.xlabel("Datetime")
plt.ylabel("CO (mg/m³)")
plt.legend()
plt.tight_layout()
plt.show()

이상치 개수: 215
상한: 5.60, 하한: -1.60

# 앞 예제의 마지막에 이어서 코딩
median = series.median()
series_corrected = series.copy()
series_corrected[outliers.index] = median

plt.figure(figsize=(12,5))
plt.plot(series.index, series, label="Before Correction", color="gray", alpha=0.6)
plt.plot(series_corrected.index, series_corrected, label="After Correction", color="blue", linewidth=1.2)
plt.scatter(outliers.index, outliers, color="red", label="Detected Outliers")
plt.title("Outlier Correction (Median Replacement)")
plt.xlabel("Datetime")
plt.ylabel("CO (mg/m³)")
plt.legend()
plt.tight_layout()
plt.show()

# ============================================
# Air Quality UCI 시계열 종합 실습 코드
# - 파일: AirQualityUCI.xlsx
# - 전처리: Date/Time → Datetime, asfreq('h'), -200→NaN, 보간+ffill/bfill
# - 분석: 장기추세, 계절별 변화, 일중패턴, 센서 상관(히트맵)
# ============================================

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# -------------------------------
# 1) 데이터 로드 + Datetime 정규화
# -------------------------------
# (필요 패키지) pip install openpyxl seaborn
df = pd.read_excel("AirQualityUCI.xlsx", sheet_name=0, engine="openpyxl")

# Time이 HH.MM.SS 형태인 경우가 있어 '.' → ':'로 교체, 공백/NaN은 00:00:00으로 대체
df["Time"] = (
    df["Time"].astype(str)
    .str.replace(".", ":", n=2, regex=False)
    .str.strip()
    .where(lambda s: s.str.len() > 0, "00:00:00")
)

# Date를 datetime으로 파싱 후 문자열화(YYYY-MM-DD) → Time과 안전하게 결합
date_ser = pd.to_datetime(df["Date"], dayfirst=True, errors="coerce")
date_str = date_ser.dt.strftime("%Y-%m-%d")
df["Datetime"] = pd.to_datetime(date_str + " " + df["Time"], errors="coerce")

# Datetime 인덱스 설정 + 시간 주기 명시 (h)
df = df.set_index("Datetime").sort_index().asfreq("h")

# -------------------------------
# 2) 센서 오류값 -200 → NaN, 숫자형 컬럼만 보간
# -------------------------------
num_cols = df.select_dtypes(include=np.number).columns
df[num_cols] = df[num_cols].replace(-200, np.nan)

# 선형보간(양방향) → 남은 NaN은 ffill/bfill로 보완
df[num_cols] = df[num_cols].infer_objects(copy=False)
df[num_cols] = df[num_cols].interpolate(method="linear", limit_direction="both")
df[num_cols] = df[num_cols].ffill().bfill()

# 분석에 자주 쓰는 주요 컬럼 후보 (실제 파일에 맞춰 자동 선택)
candidates = ["CO(GT)", "NOx(GT)", "NO2(GT)", "C6H6(GT)"]
used = [c for c in candidates if c in df.columns]
if not used:
    raise KeyError("분석용 주요 컬럼(CO(GT), NOx(GT), NO2(GT), C6H6(GT))을 찾을 수 없습니다. 파일 컬럼명을 확인하세요.")

print("[사용 컬럼]", used)

# --------------------------------
# 3) 장기 추세: 30일 이동평균(rolling)
# --------------------------------
# 시간 단위가 시간(h)이므로 30일 ≈ 30*24=720시간 윈도우
roll = df[used].rolling(window=24*30, min_periods=1).mean()

plt.figure(figsize=(12, 5))
for col in used:
    plt.plot(roll.index, roll[col], label=f"{col} (30D MA)")
plt.title("Long-term Trend (30-day Moving Average)")
plt.xlabel("Datetime")
plt.ylabel("Concentration (unit varies)")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# --------------------------------
# 4) 계절별 변화: 월별 평균 (Monthly mean)
# --------------------------------
monthly_mean = df[used].resample("M").mean()

plt.figure(figsize=(12, 5))
for col in used:
    plt.plot(monthly_mean.index, monthly_mean[col], marker="o", label=col)
plt.title("Seasonality: Monthly Mean")
plt.xlabel("Month")
plt.ylabel("Concentration (unit varies)")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# --------------------------------
# 5) 일중 패턴: 시간대별 평균 (Hour-of-day)
# --------------------------------
# 하루 0~23시 기준 평균
hourly_pattern = df[used].groupby(df.index.hour).mean()

plt.figure(figsize=(12, 5))
for col in used:
    plt.plot(hourly_pattern.index, hourly_pattern[col], marker="o", label=col)
plt.xticks(range(0, 24))
plt.title("Daily Pattern: Mean by Hour of Day")
plt.xlabel("Hour of Day (0–23)")
plt.ylabel("Concentration (unit varies)")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# --------------------------------
# 6) 센서 간 상관관계: 히트맵
# --------------------------------
# 수치형 전체 상관 또는 주요컬럼만
corr = df[used].corr(method="pearson")

plt.figure(figsize=(6, 5))
sns.heatmap(corr, annot=True, fmt=".2f", cmap="coolwarm", vmin=-1, vmax=1, square=True, cbar=True)
plt.title("Sensor/Target Correlation (Pearson)")
plt.tight_layout()
plt.show()

# --------------------------------
# 7) 보너스: 결과 요약 출력
# --------------------------------
print("\n[월별 평균 상위 5행]")
print(monthly_mean.head())
print("\n[시간대별 평균 상위 5행]")
print(hourly_pattern.head())
print("\n[상관행렬]")
print(corr)

[사용 컬럼] ['CO(GT)', 'NOx(GT)', 'NO2(GT)', 'C6H6(GT)']

[월별 평균 상위 5행]
              CO(GT)     NOx(GT)     NO2(GT)   C6H6(GT)
Datetime
2004-03-31  2.249804  140.054902  101.051961   9.935104
2004-04-30  2.138194  133.404861   91.852083  10.297856
2004-05-31  1.921841  122.704301   92.790323  10.639145
2004-06-30  1.864514  117.509028   91.445139  10.348168
2004-07-31  1.709812  123.196909  102.146505  10.595534

[시간대별 평균 상위 5행]
            CO(GT)     NOx(GT)     NO2(GT)  C6H6(GT)
Datetime
0         1.846334  198.088869  102.727612  7.852139
1         1.567172  169.398402   91.692612  6.236522
2         1.262307  136.461424   78.734202  4.704003
3         1.079787  120.127575   71.472095  3.749056
4         0.991097  104.478598   64.486141  3.312346

[상관행렬]
            CO(GT)   NOx(GT)   NO2(GT)  C6H6(GT)
CO(GT)    1.000000  0.791297  0.675975  0.816538
NOx(GT)   0.791297  1.000000  0.763096  0.655011
NO2(GT)   0.675975  0.763096  1.000000  0.570459
C6H6(GT)  0.816538  0.655011  0.570459  1.000000

# 통계 패키지 설치
pip install statsmodels

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# 예시 데이터 (월별 대기오염 농도)
date_rng = pd.date_range(start='2020-01', end='2023-01', freq='ME')
data = {
    'PM10': [80,75,70,65,90,100,95,85,80,75,70,65,
             90,100,95,85,80,75,70,65,90,100,95,85,
             80,75,70,65,90,100,95,85,80,75,70,65]
}
df = pd.DataFrame(data, index=date_rng)

# -------------------------------
# 1️⃣ 가법모형 (Additive Model)
# -------------------------------
result_add = seasonal_decompose(df['PM10'], model='additive')
fig1 = result_add.plot()
plt.show()

# -------------------------------
# 2️⃣ 승법모형 (Multiplicative Model)
# -------------------------------
result_mul = seasonal_decompose(df['PM10'], model='multiplicative')
fig2 = result_mul.plot()
plt.show()

df["diff"] = df["value"].diff()

from statsmodels.tsa.stattools import adfuller

adf_result = adfuller(df["Passengers"])
print(f"ADF Statistic: {adf_result[0]:.3f}")
print(f"p-value: {adf_result[1]:.4f}")

# 1차 차분 (Trend 제거)
df["diff1"] = df["Passengers"].diff(periods=1)

# 계절차분 (Seasonal Difference): 계절성 제거
#   12는 월별 데이터에서 1년 주기(12개월) 를 의미
df["diff_seasonal"] = df["Passengers"].diff(periods=12)

# 복합 차분: 추세 + 계절성 모두 제거
df["diff_both"] = df["Passengers"].diff(periods=1).diff(periods=12)

import numpy as np
df["log_passengers"] = np.log(df["Passengers"])

# File: 08_stationarity_differencing_practice.py
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller

# 1. 데이터 로드 (AirPassengers 데이터셋)
# 이 데이터는 1949~1960년 국제 항공 여객 수 (월별, 단위: 1000명)
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=['Month'], index_col='Month')
df.columns = ['Passengers']  # 컬럼명 통일
print(df.head())

# 2. 정상성 검정 (ADF Test)
adf_result = adfuller(df["Passengers"])
print("ADF Statistic:", round(adf_result[0], 3))
print("p-value:", round(adf_result[1], 4))

if adf_result[1] < 0.05:
    print("✅ 정상 시계열 (차분 불필요)")
else:
    print("⚠️ 비정상 시계열 (차분 필요)")

# 3. 차분 수행
df["diff1"] = df["Passengers"].diff()
df["diff_seasonal"] = df["Passengers"].diff(12)

# 4. 시각화 비교
fig, axes = plt.subplots(3, 1, figsize=(10,8), sharex=True)

axes[0].plot(df["Passengers"], label="Original", color='blue')
axes[0].set_title("Original Time Series")
axes[0].legend()

axes[1].plot(df["diff1"], label="1st Difference", color='orange')
axes[1].set_title("1st Difference (Trend Removed)")
axes[1].legend()

axes[2].plot(df["diff_seasonal"], label="Seasonal Difference", color='green')
axes[2].set_title("Seasonal Difference (Period=12)")
axes[2].legend()

plt.tight_layout()
plt.show()

            Passengers
Month
1949-01-01         112
1949-02-01         118
1949-03-01         132
1949-04-01         129
1949-05-01         121
ADF Statistic: 0.815
p-value: 0.9919
⚠️ 비정상 시계열 (차분 필요)

# File: 09_sarima_forecasting_example_v2.py
# SARIMA 모델을 이용한 항공 여객 수 예측 실습 (경고 제거 버전)
# ------------------------------------------------------------
# 이 예제는 AirPassengers 데이터셋을 사용하여
# 월별 항공 여객 수를 예측하는 SARIMA 모델을 실습하는 코드이다.
# ------------------------------------------------------------

from statsmodels.tsa.statespace.sarimax import SARIMAX
import pandas as pd
import matplotlib.pyplot as plt

# ------------------------------------------------------------
# (1) 데이터 로드 및 전처리
# ------------------------------------------------------------

# AirPassengers 데이터셋: 1949~1960년 국제 항공 여객 수 (월별)
# 데이터 출처: https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"

# CSV 파일을 불러오며, 'Month' 열을 날짜(Datetime)로 변환하고
# 이를 인덱스로 설정하여 시계열 형태로 구성
df = pd.read_csv(url, parse_dates=['Month'], index_col='Month')

# 컬럼명을 'Passengers'로 변경 (기존 컬럼명: 'Passengers'와 동일하지만 명시적으로 지정)
df.columns = ['Passengers']

# 인덱스 주기(frequency) 명시
# "MS"는 Month Start (매월 1일) 주기를 의미하며,
# statsmodels에서 주기 경고(ValueWarning)를 제거하기 위해 반드시 설정 필요
df = df.asfreq('MS')        # 또는 df.index.freq = 'MS' 도 가능

# 데이터 확인 (앞부분 5개 행)
print(f"✅ Data frequency set to: {df.index.freqstr}")
print(df.head())

# ------------------------------------------------------------
# (2) SARIMA 모델 정의
# ------------------------------------------------------------
# SARIMA(p, d, q) (P, D, Q, s)
#
# ─ 비계절(Non-seasonal) 차수: order = (p, d, q)
#   • p (AR order): 비계절 자기회귀 차수
#       - 현재 시점 y_t가 과거의 y_{t-1}, y_{t-2}, ... y_{t-p}까지 
#         몇 단계의 "값"에 의존하는지
#       - p가 크면 과거 값들을 많이 참조 → 과적합 위험 ↑, PACF 절단점으로 대략 판단
#   • d (Integration, differencing): 비계절 차분 횟수
#       - 추세(비정상성) 제거를 위해 몇 번 차분할지 (보통 0 또는 1, 드물게 2)
#       - 과도한 차분은 정보 손실/과분산 유의
#   • q (MA order): 비계절 이동평균 차수
#       - 현재 시점 y_t가 과거의 오차(잔차) e_{t-1}, e_{t-2}, ... e_{t-q}에 
#         얼마만큼 의존하는지
#       - q는 ACF 절단점으로 대략 판단
#
# ─ 계절(Seasonal) 차수: seasonal_order = (P, D, Q, s)
#   • P (Seasonal AR order): 계절 자기회귀 차수
#       - y_t가 계절 간격 s만큼 떨어진 
#         과거 값(y_{t-s}, y_{t-2s}, …, y_{t-Ps})에 의존하는 정도
#       - 예) s=12(월별 데이터)일 때, 1년 전/2년 전 같은 달의 값 의존
#   • D (Seasonal differencing): 계절 차분 횟수
#       - 계절적 평균 이동을 제거하기 위해 몇 번 계절 차분할지 (보통 0 또는 1)
#       - D=1이면 ∇_s y_t = y_t − y_{t−s} (예: 이번 달 − 작년 같은 달)
#   • Q (Seasonal MA order): 계절 이동평균 차수
#       - 현재 y_t가 계절 간격 s로 떨어진 과거 오차(e_{t-s}, e_{t-2s}, …)에 
#         의존하는 정도
#   • s (Seasonal period): 계절 주기 길이
#       - 월별 데이터의 연간 주기: s=12
#       - 분기 데이터의 연간 주기: s=4
#       - 일별 데이터의 주간 주기: s=7 등
#
# ─ 기타 옵션
#   • enforce_stationarity=False
#       - 모수 추정 시 “정상성(Stationarity) 제약”을 강제로 두지 않음
#       - 데이터가 약간 비정상이어도 추정이 수렴하도록 유연성 부여
#   • enforce_invertibility=False
#       - “가역성(Invertibility) 제약”을 강제로 두지 않음
#       - MA 성분이 경계에 가까운 경우에도 추정 허용(수렴성/안정성 트레이드오프)
#
# 예시: (1,1,1)(1,1,1,12)는
#   → 비계절 AR(1) + 비계절 차분 1회(d=1) + 비계절 MA(1)
#   → 12개월 주기의 계절 AR(1) + 계절 차분 1회(D=1) + 계절 MA(1)
#   → 즉, 추세와 12개월 계절성이 모두 존재한다고 보고 이를 동시에 모형화

non_seasonal_order = (1, 1, 1)      # (p, d, q)
seasonal_order = (1, 1, 1, 12)      # (P, D, Q, s)  ← 월별 데이터의 연간 주기

model = SARIMAX(
    endog=df["Passengers"],          # 종속변수 (시계열)
    order=non_seasonal_order,        # 비계절 차수 (p, d, q)
    seasonal_order=seasonal_order,   # 계절 차수 (P, D, Q, s)
    enforce_stationarity=False,      # 정상성 강제 X (경계 모수 허용)
    enforce_invertibility=False,     # 가역성 강제 X (경계 모수 허용)
)

# ------------------------------------------------------------
# (3) 모델 학습 (Fitting)
# ------------------------------------------------------------

# 모델 학습을 수행하면 내부적으로 MLE(Maximum Likelihood Estimation)를 통해
# AR, MA, 계절 항의 계수를 추정함
result = model.fit()

# ------------------------------------------------------------
# 결과 요약 출력 (모델 통계 요약표)
# ------------------------------------------------------------
# .summary() 출력 항목 주요 해석:
#
# [1] AIC (Akaike Information Criterion)
#     - 모델의 적합도와 복잡도(파라미터 수)를 동시에 고려하는 지표
#     - 값이 작을수록 데이터에 더 잘 맞고 과적합 위험이 낮음
#     - 여러 모델 비교 시, AIC가 가장 작은 모델 선택
#
# [2] BIC (Bayesian Information Criterion)
#     - AIC와 유사하지만, 복잡한 모델에 더 강한 패널티 부여
#     - 데이터가 많을 때 BIC가 더 보수적인 모델을 선호함
#     - 역시 값이 작을수록 좋은 모델
#
# [3] coef (Coefficient)
#     - 각 AR, MA, Seasonal, Trend 항의 추정 계수(모델 파라미터)
#     - 계수의 부호(+, -)로 변수의 영향 방향을 해석할 수 있음
#
# [4] std err (Standard Error)
#     - 해당 계수의 추정치가 얼마나 불확실한지 나타냄
#     - 작을수록 신뢰도 높은 추정 (즉, 변동성 적음)
#
# [5] z (z-statistic)
#     - 계수 유의성 검정 통계량 (계수 / 표준오차)
#     - |z| 값이 크면, 계수가 통계적으로 유의할 가능성이 높음
#
# [6] P>|z| (p-value)
#     - 유의확률: 귀무가설(H0: 계수=0)을 기각할 확률
#     - p < 0.05 → 해당 계수가 유의하게 0이 아님 (즉, 유효한 변수)
#
# [7] Ljung-Box (Q)
#     - 잔차가 백색잡음(독립적)인지 검정
#     - p > 0.05 → 잔차에 자기상관 없음 → 모델 적합 양호
#
# [8] Jarque-Bera (JB)
#     - 잔차의 정규성(평균0, 대칭성)을 검정
#     - p > 0.05 → 정규분포 가정 만족
#
# [9] Heteroskedasticity (H)
#     - 잔차의 등분산성 여부 판단
#     - 값이 1에 가까울수록 등분산(좋음)
#
# [10] Durbin–Watson
#     - 잔차의 자기상관(특히 1차 자기상관) 여부 점검
#     - 2에 가까울수록 자기상관 없음 (1 미만: 양의 상관, 3 이상: 음의 상관)
#
# 전체적으로:
#  → AIC/BIC가 작고, 잔차 진단(Ljung-Box, JB, H, DW)이 양호하면
#     모델이 데이터 패턴을 잘 설명한다고 해석할 수 있음.
#
print(result.summary())

# ------------------------------------------------------------
# (4) 미래 12개월(1년) 예측
# ------------------------------------------------------------

# 향후 12개월 데이터를 예측 (forecast)
forecast = result.get_forecast(steps=12)

# 예측 평균값
forecast_mean = forecast.predicted_mean

# 예측 구간 (95% 신뢰구간)
forecast_ci = forecast.conf_int()

# ------------------------------------------------------------
# (5) 예측 결과 시각화
# ------------------------------------------------------------

plt.figure(figsize=(10, 5))

# 실제 관측값 (1949~1960)
plt.plot(df.index, df["Passengers"], label="Observed", color='blue')

# 예측된 평균값 (1961~1962)
plt.plot(forecast_mean.index, forecast_mean, label="Forecast", color='orange')

# 예측 신뢰구간 (Confidence Interval)
plt.fill_between(
    forecast_ci.index,
    forecast_ci.iloc[:, 0],
    forecast_ci.iloc[:, 1],
    color='orange',
    alpha=0.3,
    label="Confidence Interval"
)

# 그래프 제목 및 레이블
plt.title("SARIMA Forecast for AirPassengers (1 year ahead)", fontsize=13)
plt.xlabel("Year")
plt.ylabel("Number of Passengers (in thousands)")

# 범례 및 레이아웃 조정
plt.legend()
plt.tight_layout()
plt.show()

✅ Data frequency set to: MS
            Passengers
Month
1949-01-01         112
1949-02-01         118
1949-03-01         132
1949-04-01         129
1949-05-01         121
                            SARIMAX Results
=========================================================================
Dep. Variable:                     Passengers   No. Observations: 144
Model:         SARIMAX(1, 1, 1)x(1, 1, 1, 12)   Log Likelihood   -456.103
Date:                        Fri, 17 Oct 2025   AIC               922.205
Time:                                12:24:02   BIC               936.016
Sample:                            01-01-1949   HQIC              927.812
                                 - 12-01-1960
Covariance Type:                          opg
=========================================================================
                coef    std err          z     P>|z|     [0.025    0.975]
-------------------------------------------------------------------------
ar.L1        -0.2298      0.401     -0.573     0.567     -1.016     0.557
ma.L1        -0.0987      0.374     -0.264     0.792     -0.832     0.634
ar.S.L12     -0.5460      0.299     -1.825     0.068     -1.133     0.041
ma.S.L12      0.3959      0.352      1.125     0.261     -0.294     1.086
sigma2      140.2945     17.997      7.795     0.000    105.020   175.569
=========================================================================
Ljung-Box (L1) (Q):             0.00       Jarque-Bera (JB):         5.42
Prob(Q):                        0.95       Prob(JB):                 0.07
Heteroskedasticity (H):         2.51       Skew:                     0.12
Prob(H) (two-sided):            0.01       Kurtosis:                 4.03
=========================================================================