Prognoza cen surowcow.Rmd

---
title: "Prognoza cen surowców"
author: "Agnieszka Wrzos"
output: 
  html_document:
    #toc: no
  pdf_document:
    keep_tex: yes
    #number_sections: yes
    #toc: no
  word_document:
    #toc: no
header-includes: \usepackage[polish]{babel}
                \usepackage{makeidx}
                \usepackage{amssymb}
                \usepackage{amsfonts}
                \usepackage{amsmath}
                \usepackage[T1]{fontenc}
                \usepackage{color}
                
fontsize: 12pt
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=TRUE, message=FALSE, warning=FALSE)

```

# Pallad (Palladium)

## Odczyt 

```{r, echo=FALSE}
library(tseries)
library(ggplot2)
library(dplyr)

palladium <- get.hist.quote(instrument = "PA=F", provider = "yahoo",
                      quote ="Close", start = "2015-01-01",
                      end="2021-01-30")

df_szereg <- data.frame(szereg=as.numeric(palladium), Data=time(palladium))

df_szereg <- na.omit(df_szereg)

dd <- "2021-01-20"

df_szereg1 <- df_szereg %>% filter(Data < dd)


  styl <- theme_light()+
  theme(panel.grid = element_line(colour = "#008000", linetype = "dotted"))+
  theme(panel.background = element_rect(fill = "white" , colour = "black", size=3))+
  theme(panel.border = element_rect(linetype = "solid",  size=1.2, fill = NA))+
  theme(plot.title = element_text(size=14, face="bold"))+
  theme(axis.title.y = element_text(size=11, face="bold"))+
  theme(axis.title.x = element_text(size=11, face="bold"))+
  theme(axis.text = element_text(size=9, face="bold"))+
  theme(axis.text.y = element_text(angle = 15))
```

Wczytujemy dane z Yahoo Finance od `r df_szereg$Data[1]` do `r dd`. Na podstawie tych notowań dokonano identyfikacji zachowania cen Palladium. 

## acf, pacf, wykres

```{r, echo=FALSE, fig.cap="Notowania cen palladu oraz wykresy ACF i PACF", fig.height=4, fig.width=6}
wykres_szereg <- ggplot(df_szereg1, mapping = aes(x=Data, y=szereg))+
  geom_line(size=1.05)+ggtitle("Notowania cen palladu") + 
  xlab("Data")+ylab("Cena [$]")+styl

#wykres_szereg

acf_val <- acf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
acf_df <- with(acf_val, data.frame(lag, acf))

acf_plt <- ggplot(data = acf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl
#acf_plt

pacf_val <- pacf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
pacf_df <- with(pacf_val, data.frame(lag, acf))

pacf_plt <- ggplot(data = pacf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) +
  ggtitle("PACF") + 
  xlab("Lag")+ylab("")+styl
#pacf_plt


library(grid)
library(gridExtra)

grid.arrange(wykres_szereg, acf_plt, pacf_plt,
             nrow=2, layout_matrix=rbind(c(1,1), c(2,3)))


```

## Współczynnik Hursta

```{r, echo=FALSE}
palladium <- df_szereg1$szereg
hurst <- pracma::hurstexp(palladium) # H>0.5
```

Wartość współczynnika  Hursta wynosi `r hurst$Hs`.

## Testy pierwiastka jednostkowego.


```{r}
adf.test(palladium) # czyli stopień integracji jest różny od 0

kpss.test(palladium)

pp.test(palladium)

```


## Dopasowanie modelu ARFIMA

```{r, echo=FALSE}
library(forecast)
fit1 <- arfima(palladium, drange = c(0, 0.5), estim = "mle")
tab <- summary(fit1)
```

Do oszacowania parametrów wykorzystano Metodę Najwiekszej Wiarygodności
Logarytm z funkcji wiarygodności wynosi `r tab$log.likelihood`.

```{r, echo=FALSE}

tab_cof <- as.data.frame(tab$coefficients)
tab_cof[,4] <- as.character(formatC(tab$coefficients[,4], format="e", digits = 3))


knitr::kable( tab_cof, digits = 6, caption = "Wartości parametrów dla modelu ARFIMA")
```

Wykres przedstawia wartości rzeczywiste oraz wartości dopasowane oszacowane na podstawie modelu ARFIMA

```{r, echo=FALSE, fig.cap="Dopasowanie cen Palladium za pomocą ARFIMA" , fig.height=4, fig.width=6}
df_dopasowanie <- data.frame(Data=rep(df_szereg1$Data,2),
                             y=c(df_szereg1$szereg, fit1$fitted),
                             Type=c(rep("rzeczywiste", length(df_szereg1$szereg)),
                                    rep("dopasowane", length(df_szereg1$szereg))))


ggplot(df_dopasowanie %>% filter(Data >= "2017-01-01"), mapping = aes(x=Data, y=y, colour=Type))+
  geom_line(size=1.05)+ggtitle("Notowania cen palladu") + 
  xlab("Data")+ylab("Cena [$]")+styl+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "dopasowane"),
                     values=c("#000080", "#FF0000"))

```


## Analiza reszt


```{r, echo=FALSE}
reszty <- fit1$residuals


dh <- hist(reszty, plot = FALSE, breaks = 20)
dh_data <- with(dh, data.frame(mids, density))

t <- seq(min(reszty), max(reszty), length.out = 200)
df_reszty <- data.frame(x=t, y=dnorm(t, mean = mean(reszty),
                                     sd=sd(reszty)))

r1 <- ggplot(dh_data)+geom_col(aes(x=mids, y=density), fill="grey")+
  geom_line(data=df_reszty, mapping = aes(x=x, y=y), 
            colour="red", size=1.05)+xlab("x")+ylab("Gęstość")+styl

tx <- sort(reszty)
nr <- length(reszty)
df_cdf <- data.frame(x=rep(tx,2), 
                     y=c(cumsum(rep(1/nr, nr)), 
                         pnorm(tx, mean = mean(reszty), sd=sd(reszty))),
                     type=c(rep("empiryczny", nr), rep("teoretyczny", nr)))


r2 <- ggplot(df_cdf, mapping = aes(x=x, y=y, colour=type))+
  geom_line(size=1.05)+xlab("x")+ylab("Prawdopodobieństwo")+styl+
  scale_color_manual(name="Rozkład",
                     breaks = c("empiryczny", "teoretyczny"),
                     values=c("#000080", "#FF0000"))

acf_reszty <- acf(reszty, plot = FALSE, lag.max = 50)
acf_df_res <- with(acf_reszty, data.frame(lag, acf))

acf_plt_res <- ggplot(data = acf_df_res, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl

```


```{r, echo=FALSE, fig.cap="Identyfikacja reszt", fig.height=4, fig.width=6}

grid.arrange(r1, r2, acf_plt_res, nrow=2, 
             layout_matrix = rbind(c(1,1, 2, 2, 2) ,c(3,3,3, 3,3)))
```

## Badanie normalności rozkładu reszt

```{r, echo=FALSE}
library(nortest)
test_shapiro <- shapiro.test(reszty)
test_lillie <- lillie.test(reszty)
test_ad <- ad.test(reszty)
test_cvm <- cvm.test(reszty)
Test <- c("Shapiro-Wilka", "Lillieforsa", "Andersona-Darlinga", "Cramera von Misesa")
Statystyka <- c(test_shapiro$statistic, test_lillie$statistic, test_ad$statistic, test_cvm$statistic)

p.val <- c(test_shapiro$p.value, test_lillie$p.value, test_ad$p.value, test_cvm$p.value)
p.val <- formatC(p.val, format="e", digits = 3)

test_tab <- data.frame(Test, Statystyka, p.val)

knitr::kable(test_tab, caption = "Wyniki testów badania normalności rozkładu reszt.")

```

\newpage

## Analiza niezależności reszt 

```{r, echo=FALSE}

tab_Box_test <- NULL
for (j in 1:10) {
  res_Box_test <- Box.test(reszty, lag=j, type="Ljung-Box")
  tab_Box_test <- rbind(tab_Box_test, c(j, acf_reszty$acf[j+1],
                                        res_Box_test$statistic,
                                        res_Box_test$p.value))
}

knitr::kable(tab_Box_test, col.names = c("$\\tau$", "$\\rho^{\\varepsilon}_{\\tau}$", "$\\chi^2_{(\\tau)}$", "p.val"))

```


## Predykcja

 
```{r, echo=FALSE}
library(forecast)
dl <- "2020-05-01"
predykcja <- forecast(fit1, h=dim(df_szereg)[1]-dim(df_szereg1)[1], 
                      level=c(80, 95))


df_szereg$Type <- rep("rzeczywiste", length(df_szereg$szereg))
df_szereg2 <- df_szereg %>% filter(Data >= dd)

df_predict <- data.frame(Data=df_szereg2$Data,
                         szereg=predykcja$mean)
df_predict$Type <- rep("prognozowane", length(df_predict$szereg))

df <- rbind(df_szereg, df_predict)

df_interval <- data.frame(Data=df_szereg2$Data, predykcja$lower, predykcja$upper)

names(df_interval) <- c("Data", "l80", "l95", "u80", "u95")

df_tab <- data.frame(df_interval, df_szereg2$szereg, df_predict$szereg)
names(df_tab)[6:7] <- c("rzeczywiste", "prognozowane")  
knitr::kable(df_tab)

  
```


Rysunek poniżej przedstawia wartoći notowania cen Palladium od `r dl` do `r dd` - kolor granatowy. W celu weryfikacji prognoz również na wykresie zaznaczono wartoći rzeczywiste (kolor granatowy) oraz prognozowane (kolor czerwony). Dodatkowo wynaczono przedizały ufności dla poziomu $80 \%$ oraz $95 \%$.


```{r, echo=FALSE, , fig.height=4, fig.width=6}
ggplot()+
  geom_ribbon(data = df_interval, aes(Data, ymin = l95, ymax = u95, fill = "95%")) +
  geom_ribbon(data = df_interval, aes(Data, ymin = l80, ymax = u80, fill = "80%"))+
  geom_line(data=df %>% filter(Data >= dl), 
                   mapping = aes(x=Data, y=szereg, colour=Type), size=1.05)+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "prognozowane"),
                     values=c("#000080", "#FF0000"))+
  scale_fill_manual(name = "Prognoza przedziałowa",
                    values = c("95%" = "#BDB76B", "80%" = "#FFEFD5")) +
  guides(colour = guide_legend(order = 1), fill = guide_legend(order = 2))+
  xlab("Data")+ylab("Cena [$]")+styl
```



# Kakao (Cocoa)

## Odczyt 

```{r, echo=FALSE}
library(tseries)
library(ggplot2)
library(dplyr)

cocoa <- get.hist.quote(instrument = "CC=F", provider = "yahoo",
                      quote ="Close", start = "2015-01-01",
                      end="2021-01-30")

df_szereg <- data.frame(szereg=as.numeric(cocoa), Data=time(cocoa))

df_szereg <- na.omit(df_szereg)

dd <- "2021-01-20"

df_szereg1 <- df_szereg %>% filter(Data < dd)


  styl <- theme_light()+
  theme(panel.grid = element_line(colour = "#008000", linetype = "dotted"))+
  theme(panel.background = element_rect(fill = "white" , colour = "black", size=3))+
  theme(panel.border = element_rect(linetype = "solid",  size=1.2, fill = NA))+
  theme(plot.title = element_text(size=14, face="bold"))+
  theme(axis.title.y = element_text(size=11, face="bold"))+
  theme(axis.title.x = element_text(size=11, face="bold"))+
  theme(axis.text = element_text(size=9, face="bold"))+
  theme(axis.text.y = element_text(angle = 15))
```

Wczytujemy dane z ... od `r df_szereg$Data[1]` do `r dd`. Na podstawie tych notowań dokonano identyfikacji zachowania cen kakao. 

## acf, pacf, wykres

```{r, echo=FALSE, fig.cap="Notowania cen kakao oraz wykresy ACF i PACF", fig.height=4, fig.width=6}
wykres_szereg <- ggplot(df_szereg1, mapping = aes(x=Data, y=szereg))+
  geom_line(size=1.05)+ggtitle("Notowania cen kakao") + 
  xlab("Data")+ylab("Cena [$]")+styl

#wykres_szereg

acf_val <- acf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
acf_df <- with(acf_val, data.frame(lag, acf))

acf_plt <- ggplot(data = acf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl
#acf_plt

pacf_val <- pacf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
pacf_df <- with(pacf_val, data.frame(lag, acf))

pacf_plt <- ggplot(data = pacf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) +
  ggtitle("PACF") + 
  xlab("Lag")+ylab("")+styl
#pacf_plt


library(grid)
library(gridExtra)

grid.arrange(wykres_szereg, acf_plt, pacf_plt,
             nrow=2, layout_matrix=rbind(c(1,1), c(2,3)))


```

## Współczynnik Hursta


```{r, echo=FALSE}
cocoa <- df_szereg1$szereg
hurst <- pracma::hurstexp(cocoa) # H>0.5
```

Wartość współczynnika  Hursta wynosi `r hurst$Hs`.

## Testy pierwiastka jednostkowego

```{r}
adf.test(cocoa) # czyli stopień integracji jest różny od 0

kpss.test(cocoa)

pp.test(cocoa)

```


## Dopasowanie modelu ARFIMA

```{r, echo=FALSE}
library(forecast)
fit1 <- arfima(cocoa, drange = c(0, 0.5), estim = "mle")
tab <- summary(fit1)
```

Do oszacowania parametrów wykorzystano Metodę Najwiekszej Wiarygodności
Logarytm z funkcji wiarygodności wynosi `r tab$log.likelihood`, natomiast tabela poniżej podaje wartości parametrów strukturalnych oraz wyniki testu istotności tych parametrów

```{r, echo=FALSE}

tab_cof <- as.data.frame(tab$coefficients)
tab_cof[,4] <- as.character(formatC(tab$coefficients[,4], format="e", digits = 3))


knitr::kable( tab_cof, digits = 6, caption = "Wartości parametrów dla modelu ARFIMA")
```

Wykres przedstawia wartości rzeczywiste oraz wartości dopasowane oszacowane na podstawie modelu ARFIMA

```{r, echo=FALSE, fig.cap="Dopasowanie cen kakao za pomocą ARFIMA" , fig.height=4, fig.width=6}
df_dopasowanie <- data.frame(Data=rep(df_szereg1$Data,2),
                             y=c(df_szereg1$szereg, fit1$fitted),
                             Type=c(rep("rzeczywiste", length(df_szereg1$szereg)),
                                    rep("dopasowane", length(df_szereg1$szereg))))


ggplot(df_dopasowanie %>% filter(Data >= "2017-01-01"), mapping = aes(x=Data, y=y, colour=Type))+
  geom_line(size=1.05)+ggtitle("Notowania cen kakao") + 
  xlab("Data")+ylab("Cena [$]")+styl+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "dopasowane"),
                     values=c("#000080", "#FF0000"))

```



## Analiza reszt


```{r, echo=FALSE}
reszty <- fit1$residuals


dh <- hist(reszty, plot = FALSE, breaks = 20)
dh_data <- with(dh, data.frame(mids, density))

t <- seq(min(reszty), max(reszty), length.out = 200)
df_reszty <- data.frame(x=t, y=dnorm(t, mean = mean(reszty),
                                     sd=sd(reszty)))

r1 <- ggplot(dh_data)+geom_col(aes(x=mids, y=density), fill="grey")+
  geom_line(data=df_reszty, mapping = aes(x=x, y=y), 
            colour="red", size=1.05)+xlab("x")+ylab("Gęstość")+styl

tx <- sort(reszty)
nr <- length(reszty)
df_cdf <- data.frame(x=rep(tx,2), 
                     y=c(cumsum(rep(1/nr, nr)), 
                         pnorm(tx, mean = mean(reszty), sd=sd(reszty))),
                     type=c(rep("empiryczny", nr), rep("teoretyczny", nr)))


r2 <- ggplot(df_cdf, mapping = aes(x=x, y=y, colour=type))+
  geom_line(size=1.05)+xlab("x")+ylab("Prawdopodobieństwo")+styl+
  scale_color_manual(name="Rozkład",
                     breaks = c("empiryczny", "teoretyczny"),
                     values=c("#000080", "#FF0000"))

acf_reszty <- acf(reszty, plot = FALSE, lag.max = 50)
acf_df_res <- with(acf_reszty, data.frame(lag, acf))

acf_plt_res <- ggplot(data = acf_df_res, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl

```


```{r, echo=FALSE, fig.cap="Identyfikacja reszt", fig.height=4, fig.width=6}

grid.arrange(r1, r2, acf_plt_res, nrow=2, 
             layout_matrix = rbind(c(1,1, 2, 2, 2) ,c(3,3,3, 3,3)))
```



```{r, echo=FALSE}
library(nortest)
test_shapiro <- shapiro.test(reszty)
test_lillie <- lillie.test(reszty)
test_ad <- ad.test(reszty)
test_cvm <- cvm.test(reszty)
Test <- c("Shapiro-Wilka", "Lillieforsa", "Andersona-Darlinga", "Cramera von Misesa")
Statystyka <- c(test_shapiro$statistic, test_lillie$statistic, test_ad$statistic, test_cvm$statistic)

p.val <- c(test_shapiro$p.value, test_lillie$p.value, test_ad$p.value, test_cvm$p.value)
p.val <- formatC(p.val, format="e", digits = 3)

test_tab <- data.frame(Test, Statystyka, p.val)

knitr::kable(test_tab, caption = "Wyniki testów badania normalności rozkładu reszt.")

```

\newpage

## Analiza niezależności reszt 

```{r, echo=FALSE}

tab_Box_test <- NULL
for (j in 1:10) {
  res_Box_test <- Box.test(reszty, lag=j, type="Ljung-Box")
  tab_Box_test <- rbind(tab_Box_test, c(j, acf_reszty$acf[j+1],
                                        res_Box_test$statistic,
                                        res_Box_test$p.value))
}

knitr::kable(tab_Box_test, col.names = c("$\\tau$", "$\\rho^{\\varepsilon}_{\\tau}$", "$\\chi^2_{(\\tau)}$", "p.val"))

```

## Predykcja

 
```{r, echo=FALSE}
library(forecast)
dl <- "2020-05-01"
predykcja <- forecast(fit1, h=dim(df_szereg)[1]-dim(df_szereg1)[1], 
                      level=c(80, 95))


df_szereg$Type <- rep("rzeczywiste", length(df_szereg$szereg))
df_szereg2 <- df_szereg %>% filter(Data >= dd)

df_predict <- data.frame(Data=df_szereg2$Data,
                         szereg=predykcja$mean)
df_predict$Type <- rep("prognozowane", length(df_predict$szereg))

df <- rbind(df_szereg, df_predict)

df_interval <- data.frame(Data=df_szereg2$Data, predykcja$lower, predykcja$upper)

names(df_interval) <- c("Data", "l80", "l95", "u80", "u95")


  
```


Rysunek poniżej przedstawia wartoći notowania cen kakao od `r dl` do `r dd` - kolor granatowy. W celu weryfikacji prognoz również na wykresie zaznaczono wartoći rzeczywiste (kolor granatowy) oraz prognozowane (kolor czerwony). Dodatkowo wynaczono przedizały ufności dla poziomu $80 \%$ oraz $95 \%$.


```{r, echo=FALSE, , fig.height=4, fig.width=6}
ggplot()+
  geom_ribbon(data = df_interval, aes(Data, ymin = l95, ymax = u95, fill = "95%")) +
  geom_ribbon(data = df_interval, aes(Data, ymin = l80, ymax = u80, fill = "80%"))+
  geom_line(data=df %>% filter(Data >= dl), 
                   mapping = aes(x=Data, y=szereg, colour=Type), size=1.05)+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "prognozowane"),
                     values=c("#000080", "#FF0000"))+
  scale_fill_manual(name = "Prognoza przedziałowa",
                    values = c("95%" = "#BDB76B", "80%" = "#FFEFD5")) +
  guides(colour = guide_legend(order = 1), fill = guide_legend(order = 2))+
  xlab("Data")+ylab("Cena [$]")+styl
```




# Olej opałowy (Heating Oil)

## odczyt 

```{r, echo=FALSE}
library(tseries)
library(ggplot2)
library(dplyr)

heating_oil <- get.hist.quote(instrument = "HO=F", provider = "yahoo",
                      quote ="Close", start = "2015-01-01",
                      end="2021-01-30")

df_szereg <- data.frame(szereg=as.numeric(heating_oil), Data=time(heating_oil))

df_szereg <- na.omit(df_szereg)

dd <- "2021-01-20"

df_szereg1 <- df_szereg %>% filter(Data < dd)


  styl <- theme_light()+
  theme(panel.grid = element_line(colour = "#008000", linetype = "dotted"))+
  theme(panel.background = element_rect(fill = "white" , colour = "black", size=3))+
  theme(panel.border = element_rect(linetype = "solid",  size=1.2, fill = NA))+
  theme(plot.title = element_text(size=14, face="bold"))+
  theme(axis.title.y = element_text(size=11, face="bold"))+
  theme(axis.title.x = element_text(size=11, face="bold"))+
  theme(axis.text = element_text(size=9, face="bold"))+
  theme(axis.text.y = element_text(angle = 15))
```

## acf, pacf, wykres

```{r, echo=FALSE, fig.cap="Notowania cen oleju opałowego oraz wykresy ACF i PACF", fig.height=4, fig.width=6}
wykres_szereg <- ggplot(df_szereg1, mapping = aes(x=Data, y=szereg))+
  geom_line(size=1.05)+ggtitle("Notowania cen oleju opałowego") + 
  xlab("Data")+ylab("Cena [$]")+styl

#wykres_szereg

acf_val <- acf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
acf_df <- with(acf_val, data.frame(lag, acf))

acf_plt <- ggplot(data = acf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl
#acf_plt

pacf_val <- pacf(df_szereg1$szereg, plot = FALSE, lag.max = 100)
pacf_df <- with(pacf_val, data.frame(lag, acf))

pacf_plt <- ggplot(data = pacf_df, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) +
  ggtitle("PACF") + 
  xlab("Lag")+ylab("")+styl
#pacf_plt


library(grid)
library(gridExtra)

grid.arrange(wykres_szereg, acf_plt, pacf_plt,
             nrow=2, layout_matrix=rbind(c(1,1), c(2,3)))


```

## Współczynnik Hursta

```{r, echo=FALSE}
heating_oil <- df_szereg1$szereg
hurst <- pracma::hurstexp(heating_oil) # H>0.5
```

Wartość współczynnika  Hursta wynosi `r hurst$Hs`.

## Testy pierwiastka jednostkowego 

```{r}
adf.test(heating_oil) 

kpss.test(heating_oil)

pp.test(heating_oil)

```


## Dopasowanie modelu ARFIMA

```{r, echo=FALSE}
library(forecast)
fit1 <- arfima(heating_oil, drange = c(0, 0.5), estim = "mle")
tab <- summary(fit1)
```
Do oszacowania parametrów wykorzystano Metodę Najwiekszej Wiarygodności
Logarytm z funkcji wiarygodności wynosi `r tab$log.likelihood`

## Szacowanie parametrów MNK

```{r, echo=FALSE}

tab_cof <- as.data.frame(tab$coefficients)
tab_cof[,4] <- as.character(formatC(tab$coefficients[,4], format="e", digits = 3))


knitr::kable( tab_cof, digits = 6, caption = "Wartości parametrów dla modelu ARFIMA")
```



```{r, echo=FALSE, fig.cap="Dopasowanie cen oleju opałowego za pomocą ARFIMA" , fig.height=4, fig.width=6}
df_dopasowanie <- data.frame(Data=rep(df_szereg1$Data,2),
                             y=c(df_szereg1$szereg, fit1$fitted),
                             Type=c(rep("rzeczywiste", length(df_szereg1$szereg)),
                                    rep("dopasowane", length(df_szereg1$szereg))))


ggplot(df_dopasowanie %>% filter(Data >= "2017-01-01"), mapping = aes(x=Data, y=y, colour=Type))+
  geom_line(size=1.05)+ggtitle("Notowania cen oleju opałowego") + 
  xlab("Data")+ylab("Cena [$]")+styl+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "dopasowane"),
                     values=c("#000080", "#FF0000"))

```



## Analiza reszt 


```{r, echo=FALSE}
reszty <- fit1$residuals


dh <- hist(reszty, plot = FALSE, breaks = 20)
dh_data <- with(dh, data.frame(mids, density))

t <- seq(min(reszty), max(reszty), length.out = 200)
df_reszty <- data.frame(x=t, y=dnorm(t, mean = mean(reszty),
                                     sd=sd(reszty)))

r1 <- ggplot(dh_data)+geom_col(aes(x=mids, y=density), fill="grey")+
  geom_line(data=df_reszty, mapping = aes(x=x, y=y), 
            colour="red", size=1.05)+xlab("x")+ylab("Gęstość")+styl

tx <- sort(reszty)
nr <- length(reszty)
df_cdf <- data.frame(x=rep(tx,2), 
                     y=c(cumsum(rep(1/nr, nr)), 
                         pnorm(tx, mean = mean(reszty), sd=sd(reszty))),
                     type=c(rep("empiryczny", nr), rep("teoretyczny", nr)))


r2 <- ggplot(df_cdf, mapping = aes(x=x, y=y, colour=type))+
  geom_line(size=1.05)+xlab("x")+ylab("Prawdopodobieństwo")+styl+
  scale_color_manual(name="Rozkład",
                     breaks = c("empiryczny", "teoretyczny"),
                     values=c("#000080", "#FF0000"))

acf_reszty <- acf(reszty, plot = FALSE, lag.max = 50)
acf_df_res <- with(acf_reszty, data.frame(lag, acf))

acf_plt_res <- ggplot(data = acf_df_res, mapping = aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(mapping = aes(xend = lag, yend = 0)) + 
  ggtitle("ACF") + 
  xlab("Lag")+ylab("")+styl

```


```{r, echo=FALSE, fig.cap="Identyfikacja reszt", fig.height=4, fig.width=6}

grid.arrange(r1, r2, acf_plt_res, nrow=2, 
             layout_matrix = rbind(c(1,1, 2, 2, 2) ,c(3,3,3, 3,3)))
```

## Normalność rozkładu reszt

```{r, echo=FALSE}
library(nortest)
test_shapiro <- shapiro.test(reszty)
test_lillie <- lillie.test(reszty)
test_ad <- ad.test(reszty)
test_cvm <- cvm.test(reszty)
Test <- c("Shapiro-Wilka", "Lillieforsa", "Andersona-Darlinga", "Cramera von Misesa")
Statystyka <- c(test_shapiro$statistic, test_lillie$statistic, test_ad$statistic, test_cvm$statistic)

p.val <- c(test_shapiro$p.value, test_lillie$p.value, test_ad$p.value, test_cvm$p.value)
p.val <- formatC(p.val, format="e", digits = 3)

test_tab <- data.frame(Test, Statystyka, p.val)

knitr::kable(test_tab, caption = "Wyniki testów badania normalności rozkładu reszt.")

```



## Analiza niezależności reszt 

```{r, echo=FALSE}

tab_Box_test <- NULL
for (j in 1:10) {
  res_Box_test <- Box.test(reszty, lag=j, type="Ljung-Box")
  tab_Box_test <- rbind(tab_Box_test, c(j, acf_reszty$acf[j+1],
                                        res_Box_test$statistic,
                                        res_Box_test$p.value))
}

knitr::kable(tab_Box_test, col.names = c("$\\tau$", "$\\rho^{\\varepsilon}_{\\tau}$", "$\\chi^2_{(\\tau)}$", "p.val"))

```




## Predykcja

 
```{r, echo=FALSE}
library(forecast)
dl <- "2020-05-01"
predykcja <- forecast(fit1, h=dim(df_szereg)[1]-dim(df_szereg1)[1], 
                      level=c(80, 95))


df_szereg$Type <- rep("rzeczywiste", length(df_szereg$szereg))
df_szereg2 <- df_szereg %>% filter(Data >= dd)



df_predict <- data.frame(Data=df_szereg2$Data,
                         szereg=predykcja$mean)
df_predict$Type <- rep("prognozowane", length(df_predict$szereg))

df <- rbind(df_szereg, df_predict)

df_interval <- data.frame(Data=df_szereg2$Data, predykcja$lower, predykcja$upper)

names(df_interval) <- c("Data", "l80", "l95", "u80", "u95")

df_tab <- data.frame(df_interval, df_szereg2$szereg, df_predict$szereg)
names(df_tab)[6:7] <- c("rzeczywiste", "prognozowane")  
knitr::kable(df_tab)

  
```


## Wartości notowania cen oleju opałowego

```{r, echo=FALSE, , fig.height=4, fig.width=6}
ggplot()+
  geom_ribbon(data = df_interval, aes(Data, ymin = l95, ymax = u95, fill = "95%")) +
  geom_ribbon(data = df_interval, aes(Data, ymin = l80, ymax = u80, fill = "80%"))+
  geom_line(data=df %>% filter(Data >= dl), 
                   mapping = aes(x=Data, y=szereg, colour=Type), size=1.05)+
  scale_color_manual(name="Wartości",
                     breaks = c("rzeczywiste", "prognozowane"),
                     values=c("#000080", "#FF0000"))+
  scale_fill_manual(name = "Prognoza przedziałowa",
                    values = c("95%" = "#BDB76B", "80%" = "#FFEFD5")) +
  guides(colour = guide_legend(order = 1), fill = guide_legend(order = 2))+
  xlab("Data")+ylab("Cena [$]")+styl
```