#####################################################################
#
# R package kDGLM
#
# Bayesian Analysis of Dynamic Generalized Linear Models 
# 
# Provide routines for filtering and smoothing, 
# forecasting, sampling and Bayesian analysis of Dynamic 
# Generalized Linear Models.
#
# References:
#
# 1) An Efficient Sequential Approach for k-Parametric 
#    Dynamic Generalised Linear Models by Mariane Alves, 
#    Helio Migon, Silvaneo Santos Jr and Raíra Marotta.
#    https://arxiv.org/abs/2201.05387
#
# 2) kDGLM: a R package for Bayesian analysis of Generalized 
#    Dynamic Linear Models by Silvaneo Santos Jr., 
#    Mariane Alves and Helio Migon.
#    https://arxiv.org/abs/2403.13069
#
# Tutorials: 
#
# https://silvaneojunior.github.io/kDGLM/
# https://cran.r-project.org/web/packages/kDGLM/refman/kDGLM.html
#
#####################################################################

#remotes::install_github('silvaneojunior/kDGLM')
library("kDGLM")

# Poisson case
# Monthly Airline Passenger Numbers 1949-1960
data = c(AirPassengers)
par(mfrow=c(1,1),mar = c(2, 2, 1, 1))
ts.plot(data,xlab="Monthly Airline Passenger Numbers 1949-1960")

level = polynomial_block(rate = 1, order = 2, D = 0.95)
season = harmonic_block(rate = 1, order = 2, period = 12, D = 0.975)
outcome = Poisson(lambda = "rate", data = data)
fitted.data = fit_model(level,season,AirPassengers = outcome)
summary(fitted.data)
plot(fitted.data, plot.pkg = "ggplot2")

# Gamma case 
# Corn and wheat prices from 1986 to 2014
par(mfrow=c(2,2),mar = c(2, 2, 1, 1))
ts.plot(cornWheat$corn.price,ylab="",main="Corn price")
ts.plot(cornWheat$wheat.price,ylab="",main="Wheat price")
ts.plot(cornWheat$corn.log.return,ylab="",main="Corn log-return")
ts.plot(cornWheat$wheat.log.return,ylab="",main="Wheat log-return")

Y = cornWheat$corn.log.return[1:1000]
Y = abs(Y-mean(Y))
par(mfrow=c(1,1))
plot(Y,type="h")

structure = polynomial_block(mu = 1, D = 0.95) 
outcome = Gamma(phi = 0.5, mu = "mu", data = Y) 
fitted.data = fit_model(structure, corn = outcome) 
summary(fitted.data) 
plot(fitted.data, plot.pkg = "ggplot2")