#!/usr/bin/env Rscript

df <- read.table("jj.dat", sep = "\n")


df_ts <- ts(df, start = 1960, frequency = 4)
plot(df_ts, main = "df_ts", ylab = "value")

ts_diff <- diff(df_ts)
plot(ts_diff)

decomp_add <- decompose(df_ts, type="additive")
plot(decomp_add)

df_ts_log <- log10(df_ts)
plot(df_ts_log)

ts_diff_log10 <- diff(df_ts_log)
plot(ts_diff_log10)

acf(ts_diff_log10)
pacf(ts_diff_log10)


model1 <- arima(ts_diff_log10, order = c(1, 1, 0))
model2 <- arima(ts_diff_log10, order = c(0, 1, 1))
model3 <- arima(ts_diff_log10, order = c(1, 1, 1))

# 2k - 2 ln(lambda)
AIC(model1, model2, model3)


sim_ar4 <- arima.sim(n = 10000, list(ar = c(0.9, -0.5, 0.2, -0.3)))
acf(sim_ar4, main = "ACF AR(4)")
pacf(sim_ar4, main = "PACF AR(4)")
ar4 <- arima(sim_ar4, order = c(4, 0, 0))
AIC(ar4)

sim_ar3 <- arima.sim(n = 10000, list(ar = c(0.9, -0.5, 0.2)))
acf(sim_ar3, main = "ACF AR(3)")
pacf(sim_ar3, main = "PACF AR(3)")
ar3 <- arima(sim_ar3, order = c(3, 0, 0))
AIC(ar3)

sim_ar2 <- arima.sim(n = 10000, list(ar = c(0.9, -0.5)))
acf(sim_ar2, main = "ACF AR(2)")
pacf(sim_ar2, main = "PACF AR(2)")
ar2 <- arima(sim_ar2, order = c(2, 0, 0))
AIC(ar2)

sim_ar1 <- arima.sim(n = 10000, list(ar = c(0.9)))
acf(sim_ar1, main = "ACF AR(1)")
pacf(sim_ar1, main = "PACF AR(1)")
ar1 <- arima(sim_ar1, order = c(1, 0, 0))
AIC(ar1)

sim_ma4 <- arima.sim(n = 10000, list(ma = c(-1.9, 1.7, -1.5, 1.5)))
acf(sim_ma4, main = "ACF MA(4)")
pacf(sim_ma4, main = "PACF MA(4)")
ma4 <- arima(sim_ma4, order = c(4, 0, 0))
AIC(ma4)

sim_ma3 <- arima.sim(n = 10000, list(ma = c(-1.9, 1.7, -1.5)))
acf(sim_ma3, main = "ACF MA(3)")
pacf(sim_ma3, main = "PACF MA(3)")
ma3 <- arima(sim_ma3, order = c(3, 0, 0))
AIC(ma3)

sim_ma2 <- arima.sim(n = 10000, list(ma = c(-1.9, 1.7)))
acf(sim_ma2, main = "ACF MA(2)")
pacf(sim_ma2, main = "PACF MA(2)")
ma2 <- arima(sim_ma2, order = c(2, 0, 0))
AIC(ma2)

sim_ma1 <- arima.sim(n = 10000, list(ma = c(-1.9)))
acf(sim_ma1, main = "ACF MA(1)")
pacf(sim_ma1, main = "PACF MA(1)")
ma1 <- arima(sim_ma1, order = c(1, 0, 0))
AIC(ma1)