import numpy as np, pandas as pd
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

np.random.seed(123)

# AR(1) estacionário com phi=0.5
ar_params = [1, -0.5]  # [1, -phi1, -phi2, ...]
arma_p = ArmaProcess(ar=ar_params, ma=[1])
diff_y = arma_p.generate_sample(nsample=500)

# Integrar (I=1) → ARIMA(1,1,0)
y  = np.cumsum(diff_y)   # serie nao estacionaria
dy = np.diff(y)           # primeira diferenca → AR(1)

np.random.seed(123)

# MA(1) com theta=0.9
arma_p2 = ArmaProcess(ar=[1], ma=[1, 0.9])
diff_ima = arma_p2.generate_sample(nsample=500)

# Integrar → IMA(1,1)
dates = pd.date_range('2020-01-01', periods=500, freq='D')
ima1  = pd.Series(np.cumsum(diff_ima), index=dates)

from statsmodels.tsa.arima.model import ARIMA

model   = ARIMA(endog=ima1, order=(0, 1, 1))
results = model.fit()
print(results.summary())

from statsmodels.stats.diagnostic import acorr_ljungbox

res = results.resid

# Visualizar ACF dos resíduos
fig, ax = plt.subplots(figsize=(10, 4))
plot_acf(res, lags=24, zero=False, ax=ax,
         title='ACF dos Resíduos: ARIMA(0,1,1)')
plt.show()

# Teste Ljung-Box (H₀: resíduos são ruído branco)
lb = acorr_ljungbox(res, lags=[24], model_df=2, return_df=True)
# p-value > 0.05 → não rejeitamos H₀ → modelo OK ✓

# Previsão 30 passos à frente com IC 95%
forecast = results.get_forecast(steps=30)
fc_mean  = forecast.predicted_mean
fc_ci    = forecast.conf_int(alpha=0.05)

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(ima1.iloc[-80:], color='#E4002B', label='Histórico')
ax.plot(fc_mean, color='#1a9b57', linestyle='--', label='Previsão')
ax.fill_between(fc_ci.index,
                fc_ci.iloc[:,0], fc_ci.iloc[:,1],
                alpha=0.18, color='#1a9b57', label='IC 95%')
ax.axvline(ima1.index[-1], color='gray', lw=1, ls=':')
ax.legend(); plt.tight_layout(); plt.show()

from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA

sf = StatsForecast(
    models=[AutoARIMA(season_length=1)],  # sem sazonalidade
    freq='D', n_jobs=1,
)

df = pd.DataFrame({
    'unique_id': 'ima1',
    'ds': ima1.index,
    'y':  ima1.values
})

sf.fit(df=df)
previsao = sf.predict(h=30, level=[80, 95])

# Inspecionar o modelo selecionado
modelo = sf.fitted_[0, 0].model_
arma   = modelo['arma']        # [p, q, P, Q, s, d, D]
p, q, P, Q, s, d, D = arma
print(f"Modelo selecionado: ARIMA({p},{d},{q})({P},{D},{Q})[{s}]")
print(f"AIC  : {modelo['aic']:.3f}")
print(f"BIC  : {modelo['bic']:.3f}")
print(f"sigma: {modelo['sigma2']**0.5:.4f}")

from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.api as sm

data = sm.datasets.get_rdataset('AirPassengers', 'datasets').data
AP   = pd.Series(data['value'].values,
                 index=pd.date_range('1949-01', periods=144, freq='MS'))

decomp = seasonal_decompose(AP, model='multiplicative', period=12)
decomp.plot(); plt.show()

# Transformação log converte multiplicativo em aditivo
AP_log = np.log(AP)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Aditivo: amplitude constante (inadequado para AirPassengers)
# Multiplicativo/log: amplitude proporcional à tendência (correto)

AP_log  = np.log(AP[:'1959'])      # excluir 1960 para teste
AP_diff      = AP_log.diff().dropna()    # d=1
AP_diff_seas = AP_diff.diff(12).dropna() # D=1

fig, axes = plt.subplots(2, 2, figsize=(12, 6))
plot_acf( AP_diff,      ax=axes[0,0], lags=36, zero=False)
plot_pacf(AP_diff,      ax=axes[0,1], lags=36, zero=False)
plot_acf( AP_diff_seas, ax=axes[1,0], lags=36, zero=False)
plot_pacf(AP_diff_seas, ax=axes[1,1], lags=36, zero=False)
plt.show()

from itertools import product

resultados = []
for p, q in product([0, 1, 2], repeat=2):
    try:
        m = ARIMA(AP_log, order=(p,1,q),
                  seasonal_order=(0,1,0,12)).fit()
        resultados.append({'Modelo': f'SARIMA({p},1,{q})(0,1,0)',
                           'AIC': m.aic, 'AICc': m.aicc, 'BIC': m.bic})
    except: pass

pd.DataFrame(resultados).sort_values('AICc').head(4)

fit = ARIMA(AP_log[:'1959'],
            order=(1, 1, 0),
            seasonal_order=(0, 1, 0, 12)).fit()

# Previsão em escala original (exp do log)
fc_log  = fit.get_forecast(12).predicted_mean
fc_orig = np.exp(fc_log)
ci_log  = fit.get_forecast(12).conf_int()
ci_orig = np.exp(ci_log)

print(fit.summary())
print(f"\nMAE (1960): {np.abs(fc_orig - AP['1960':]).mean():.1f}")

from prophet import Prophet  # pip install prophet

df_prophet = pd.DataFrame({
    'ds': AP.index,
    'y':  AP.values
})

# Sazonalidade multiplicativa (amplitude cresce com tendência)
m = Prophet(seasonality_mode='multiplicative')
m.fit(df_prophet[df_prophet['ds'] < '1960-01-01'])

# Previsão
future   = m.make_future_dataframe(periods=12, freq='MS')
forecast = m.predict(future)

# Visualizar componentes
fig = m.plot_components(forecast)
plt.show()

# Extrair previsão para 1960 e calcular MAE
real       = AP['1960':]
fc_prophet = forecast[forecast['ds'] >= '1960-01-01']['yhat'].values
print(f"MAE Prophet: {np.abs(fc_prophet - real.values).mean():.1f}")

# Criar dummies mensais (janeiro = referência)
meses = pd.get_dummies(AP.index.month, prefix='mes', drop_first=True)
meses.index = AP.index
trend = pd.Series(range(1, len(AP)+1), index=AP.index, name='trend')
X = pd.concat([meses, trend], axis=1).astype(float)

X_train = X[:'1959']; X_test = X['1960':]

fit_arimax = ARIMA(AP[:'1959'],
                   order=(1, 1, 0),
                   exog=X_train).fit()

fc_arimax = fit_arimax.get_forecast(12, exog=X_test).predicted_mean
print(f"MAE ARIMAX: {np.abs(fc_arimax - AP['1960':]).mean():.1f}")

from prophet import Prophet
import pandas as pd

# Criar dummies de mês
df_prophet['mes'] = df_prophet['ds'].dt.month
for m in range(2, 13):  # janeiro = referência
    df_prophet[f'mes_{m}'] = (df_prophet['mes'] == m).astype(int)

regressores = [f'mes_{m}' for m in range(2, 13)]

# Modelo Prophet com dummies como regressores
m = Prophet(seasonality_mode='multiplicative')
for reg in regressores:
    m.add_regressor(reg)

m.fit(df_prophet[df_prophet['ds'] < '1960-01-01'])

# Previsão: dummies devem estar no dataframe futuro também
future = m.make_future_dataframe(periods=12, freq='MS')
future['mes'] = future['ds'].dt.month
for reg in regressores:
    future[reg] = (future['mes'] == int(reg.split('_')[1])).astype(int)

forecast = m.predict(future)

# Extrair previsão para 1960
fc_prophet_reg = forecast[forecast['ds'] >= '1960-01-01']['yhat'].values
print(f"MAE Prophet+Dummies: {np.abs(fc_prophet_reg - AP['1960':].values).mean():.1f}")

# fc_orig        → SARIMA           (slide 9)
# fc_prophet     → Prophet mult.    (slide 10)
# fc_prophet_reg → Prophet+Dummies  (slide 11)
# fc_arimax      → ARIMAX           (slide 11)

real = AP_test.values   # 12 valores reais de 1960

resultados = {
    'SARIMA':            fc_orig,
    'Prophet (mult.)':   fc_prophet,
    'Prophet + Dummies': fc_prophet_reg,
    'ARIMAX':            fc_arimax,
}

print(f"{'Modelo':25s}  {'MAE':>8}  {'MAPE':>8}")
print('-' * 48)
for modelo, fc in sorted(resultados.items(),
                          key=lambda x: np.abs(x[1] - real).mean()):
    mae_val  = np.abs(fc - real).mean()
    mape_val = np.abs((fc - real) / real * 100).mean()
    print(f"{modelo:25s}  {mae_val:8.1f}  {mape_val:7.1f}%")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
cores = ['#1a9b57', '#1a6fc4', '#cc44ff', '#cc8800']

axes[0].plot(AP.index, AP.values, color='#aaaaaa', lw=1.5, label='Histórico')
axes[0].plot(AP_test.index, real, color='#1C1C1C', lw=2.5, ls=':', label='Real 1960', zorder=6)
for (nome, fc), cor in zip(resultados.items(), cores):
    axes[0].plot(AP_test.index, fc, color=cor, lw=2, ls='--', label=nome)
axes[0].axvline(AP_train.index[-1], color='#888', lw=1, ls=':', alpha=0.7)
axes[0].set_title('Previsões: 1960'); axes[0].set_ylabel('Passageiros')
axes[0].legend(fontsize=8, loc='upper left')
plt.tight_layout(); plt.show()

vals_mae  = [np.abs(fc - real).mean()            for fc in resultados.values()]
vals_mape = [np.abs((fc - real)/real*100).mean() for fc in resultados.values()]
x = np.arange(len(resultados)); w = 0.35
fig, ax = plt.subplots(figsize=(9, 4))
ax.bar(x - w/2, vals_mae,  width=w, color=cores, alpha=0.85, label='MAE')
ax.bar(x + w/2, vals_mape, width=w, color=cores, alpha=0.5,
       edgecolor=cores, linewidth=1.5, fill=False, label='MAPE (%)')
ax.set_xticks(x); ax.set_xticklabels(resultados.keys(), fontsize=9)
ax.set_title('MAE e MAPE% por modelo'); ax.legend(); plt.tight_layout(); plt.show()