##############################################################################
#####  R skript k diplomove praci: Testy nezavislosti dvou casovych rad  #####
#####  Autor: Pavel Zdenek                                               #####
##############################################################################

#########################################
#####  Nacteni potrebnych knihoven  #####
#########################################

library(MASS)     # generovani vicerozmerneho normalniho rozdeleni
library(ggplot2)  # balicek na grafy
library(forecast) # balicek na casove rady
library(psych)    # funkce tr - trace
library(dplyr)    # funkce inner_join - prunik

#################################
#####  Pocatecni nastaveni  #####
#################################

rm(list=ls())
graphics.off()
seed = 971012
alpha = 0.05

#############################
#####  Definice funkci  #####
#############################

generate_ARMA <- function(N, m, n, AR, MA, e, Q){ # vygeneruje N + Q pozorovani ARMA(m, n) na zaklade koeficientu AR a MA, bileho sumu e, prvnich Q pozorovani pote zahodi
  Series = rep(0, (N+Q))
  for(i in 1:(N+Q)){
    Series[i] = e[i]
    if(n > 0){
      for(j in 1:n){
        if(i-j<1) break()
        Series[i] = Series[i] + MA[j]*e[i-j]
      }
    }
    if(m>0){
      for(j in 1:m){
        if(i-j<1) break()
        Series[i] = Series[i] - AR[j]*Series[i-j]
      }
    }
  }
  if(Q!=0) Series = Series[-c(1:Q)]
  return(Series)
}

cross_covariance <- function(x, y, k, WN = FALSE){ # spocita vyberovou krizovou kovarianci mezi sdruzene slabe stacionarnimi x a y v bode k. Pri WN = FALSE se ve vypoctu pouzije prumer, jinak se nepouzije.
  N = length(x)
  mean_x = mean(x)
  mean_y = mean(y)
  if(k>=0){
    x = x[1:(N - k)]
    y = y[(1 + k):N]
  }else{
    x = x[(-k + 1):N]
    y = y[1:(N + k)]
  }
  if(WN){ 
    covariance = sum(x * y)
  }
  else{ 
    covariance = sum((x - mean_x) * (y - mean_y))
  }
  covariance = covariance/N
  return(covariance)
}

cross_correlation <- function(x, y, k, WN = FALSE){ # spocita vyberovou krizovou korelaci mezi sdruzene slabe stacionarnimi x a y v bode k. Pri WN = FALSE se ve vypoctu pouzije prumer, jinak se nepouzije.
  correlation = cross_covariance(x, y, k, WN)/sqrt(cross_covariance(x, x, 0, WN)*cross_covariance(y, y, 0, WN))
  return(correlation)
}

auto_covariance <- function(x, k){ # spocita hodnotu autokovariancni funkce pro radu x v bode k
  N = length(x)
  if(k>=0){
    x1 = x[1:(N - k)]
    x2 = x[(1 + k):N]
  }else{
    x1 = x[(-k + 1):N]
    x2 = x[1:(N + k)]
  }
  covariance =sum((x1 - mean(x)) * (x2 - mean(x)))
  covariance = covariance/N
  return(covariance)
}

auto_correlation <- function(x, k){ # spocita hodnotu vyberove autokorelacni funkce pro radu x v bode k
  correlation = auto_covariance(x, k)/auto_covariance(x, 0)
  return(correlation)
}

auto_covariance_matrix <- function(x){ # spocita vyberovou autokovarincni matici pro radu x
  covariance_vector = c(1:length(x))
  for(i in 1:(length(x))){
    covariance_vector[i] = auto_covariance(x, i - 1)
  }
  covariance_vector = c(rev(covariance_vector[-1]), covariance_vector)
  covariance_matrix = covariance_vector[length(x):length(covariance_vector)]
  for(i in 1:(length(x) - 1)){
    covariance_matrix = rbind(covariance_matrix, covariance_vector[(length(x) - i):(length(covariance_vector) - i)])
  }
  return(covariance_matrix)
}

dist_cov <- function(x, y){ # spocte hodnotu vyberove distancni kovariance mezi radami x a y
  x_matrix = matrix(rep(x, length(x)), nrow = length(x) )
  distance_x = abs(x_matrix - t(x_matrix))                                                       # a_{kl}
  distance_x_sums = colSums(distance_x)/length(x)                                   
  distance_x_sums_matrix = matrix(rep(distance_x_sums, length(x)), nrow = length(x))             # a_{k.} a a_{.l}
  distance_x_average = sum(distance_x)/length(x)^2                                               # a_{..}
  final_x = distance_x - distance_x_sums_matrix - t(distance_x_sums_matrix) + distance_x_average # A_{kl}
  
  y_matrix = matrix(rep(y, length(y)), nrow = length(y) )
  distance_y = abs(y_matrix - t(y_matrix))                                                       # b_{kl}
  distance_y_sums = colSums(distance_y)/length(y)                                         
  distance_y_sums_matrix = matrix(rep(distance_y_sums, length(y)), nrow = length(y))             # b_{k.} a b_{.l}
  distance_y_average = sum(distance_y)/length(y)^2                                               # b_{..}
  final_y = distance_y - distance_y_sums_matrix - t(distance_y_sums_matrix) + distance_y_average # B_{kl}
  
  final_x_y = final_x * final_y
  dist_cov_x_y = sqrt(sum(final_x_y)/length(x)^2) 
  return(dist_cov_x_y)
}

dist_corr <- function(x, y){ # spocte hodnotu vyberove distancni korelace mezi radami x a y
  dist_corr_x_y = dist_cov(x, y)/(sqrt(dist_cov(x, x)) * sqrt(dist_cov(y, y)))
  return(dist_corr_x_y)
}

dist_test_check <- function(x, y){ # spocte testovou statistiku pro test distancni kovarianci pro rady x a y
  x_matrix = matrix(rep(x, length(x)), nrow = length(x) )
  distance_x = abs(x_matrix - t(x_matrix))                                                        # a_{kl}
  distance_x_sums = colSums(distance_x)/length(x)   
  distance_x_sums_matrix = matrix(rep(distance_x_sums, length(x)), nrow = length(x))              # a_{k.} a a_{.l}
  distance_x_average = sum(distance_x)/length(x)^2                                                # a_{..}
  final_x = distance_x - distance_x_sums_matrix - t(distance_x_sums_matrix) + distance_x_average  # A_{kl}
  
  y_matrix = matrix(rep(y, length(y)), nrow = length(y) )                                         # b_{kl}
  distance_y = abs(y_matrix - t(y_matrix))
  distance_y_sums = colSums(distance_y)/length(y)
  distance_y_sums_matrix = matrix(rep(distance_y_sums, length(y)), nrow = length(y))              # b_{k.} a b_{.l}
  distance_y_average = sum(distance_y)/length(y)^2                                                # b_{..}
  final_y = distance_y - distance_y_sums_matrix - t(distance_y_sums_matrix) + distance_y_average  # B_{kl}
  
  final_x_y = final_x * final_y
  dist_cov_x_y = sqrt(sum(final_x_y)/length(x)^2)
  
  S = sum(distance_x) * sum(distance_y)/length(x)^4         # S_N
  test = dist_cov_x_y * sqrt(length(x)/S)
  return(test)
}

compute_N_hat <- function(x, y){ # spocte N_hat pro rady x a y pro pouziti modifikovaneho t-testu nezavislosti
  covariance_matrix_x = auto_covariance_matrix(x)
  covariance_matrix_y = auto_covariance_matrix(y)
  
  B = diag(length(x))- matrix(1, nrow = length(x), ncol = length(x))/length(x)
  
  sigma_t = tr(B %*% covariance_matrix_x %*% B %*% covariance_matrix_y)/(tr(B %*% covariance_matrix_x)*tr(B %*% covariance_matrix_y))
  
  N_hat =(1/sigma_t) + 1
  return(N_hat)
}

t_ind_test_p <- function(x, y){ # vraci p-hodnotu modifikovaneho t-testu pouziteho pro rady x a y
  N_hat = compute_N_hat(x, y)
  test_statistic_t = cross_correlation(x, y, 0) * sqrt((N_hat-2)/(1-cross_correlation(x, y, 0)^2))
  p = 2 * (1 - pt(abs(test_statistic_t), N_hat-2))
  return(p)
}

dist_ind_test_p <- function(x, y){ # vraci p-hodnotu testu distancni kovarianci pro rady x a y
  p = 2 * (1 - pnorm(dist_test_check(x, y)))
  return(p)
}

Haugh_ind_test_p <- function(x, y, M, better_approx = FALSE){ # vraci p-hodnotu Haughova testu s parametrem M pro rady x a y. Pri better_approx = FALSE se pouzije testova statistika S_M, jinak S*_M
  test_statistic_Haugh = 0
  if(better_approx){
    for(i in -M:M){
      test_statistic_Haugh = test_statistic_Haugh + cross_correlation(x, y, i, WN = TRUE)^2 * length(x)  
    }
  }else{
    for(i in -M:M){
      test_statistic_Haugh = test_statistic_Haugh + length(x)^2 * cross_correlation(x, y, i, WN = TRUE)^2/(length(x) - abs(i))  
    }
  }
  p = 1 - pchisq(test_statistic_Haugh, 2*M + 1)
  return(p)
}

find_ARMA_model <- function(x, M, N){ # najde vhodny ARMA model pro radu x podle AIC. Zkousi se vsechny modely az do radu ARMA(M, N)
  min = 1/0
  for(i in 0:M){
    for(j in 0:N){
      tryCatch({
      if (i + j != 0){
        model = arima(x, order = c(i, 0, j), include.mean = FALSE)
        if(min > model$aic){
          min = model$aic
          modelx = model
        }
      }
      },error = function(e){})
    }
  }
  return(modelx)
}

find_ARMA_model_with_mean <- function(x, M, N){ # najde vhodny ARMA model s nenulovou stredni hodnotou pro radu x podle AIC. Zkousi se vsechny modely az do radu ARMA(M, N)
  min = 1/0
  for(i in 0:M){
    for(j in 0:N){
      tryCatch({
        if (i + j !=0){
          model = arima(x, order = c(i,0,j), include.mean = TRUE)
          if(min > model$aic){
            min = model$aic
            modelx = model
          }
        }
      },error = function(e){})
    }
  }
  return(modelx)
}

Simulation <- function(Number_of_tests, alpha, N, Q, mu, s1, s2, rho, M, mx, nx, ax, bx, my, ny, ay, by){ # funkce pouzita pro opakovane generovani rad a nasledne otestovani nezavislosti
  sigma <- matrix(c(s1^2, s1 * s2 * rho, s1 * s2 * rho, s2^2), 2)
  
  ### Inicializace promennych ###
  
  haugh_test_v1_rejected = 0
  haugh_test_v2_rejected = 0
  t_test_ARMA_rejected = 0
  t_test_WN_rejected = 0
  standard_t_test_ARMA_rejected = 0
  standard_t_test_WN_rejected = 0
  dist_test_rejected = 0
  model_estimation_x_right = 0
  model_estimation_y_right = 0
  
  ### Opakovane testovani ###
  
  for(i in 1:Number_of_tests){
    if(i %% 100 == 0) cat("Simulaion number:",i, "\n")
    
    ### Generovani rad x a y ###
    
    e <- mvrnorm((N+Q), mu, sigma)                          
    colnames(e) <- c("ex", "ey")
    ex <- e[ , c("ex")]                             
    ey <- e[ , c("ey")]                                   
    
    x = generate_ARMA(N, mx, nx, ax, bx, ex, Q)
    y = generate_ARMA(N, my, ny, ay, by, ey, Q)
    
    ### odhad modelů ARMA ###
    
    modelx = find_ARMA_model(x, 3, 3)
    
    modely = find_ARMA_model(y, 3, 3)
    
    if(modelx$arma[1]==mx & modelx$arma[2]==nx) model_estimation_x_right = model_estimation_x_right + 1
    
    if(modely$arma[1]==my & modely$arma[2]==ny) model_estimation_y_right = model_estimation_y_right + 1
    
    predex = modelx$residuals[1:N]
    
    predey = modely$residuals[1:N]
    
    ### Haugh test - v1 ###
    
    if(Haugh_ind_test_p(predex, predey, M, better_approx = FALSE) < alpha) haugh_test_v1_rejected = haugh_test_v1_rejected + 1
    
    ### Haugh test - v2 ###
    
    if(Haugh_ind_test_p(predex, predey, M, better_approx = TRUE) < alpha) haugh_test_v2_rejected = haugh_test_v2_rejected + 1
    
    ### modifikovany t-test - ARMA ###
    
    if(t_ind_test_p(x, y) < alpha) t_test_ARMA_rejected = t_test_ARMA_rejected + 1
    
    ### modifikovany t-test - bily sum ###
    
    if(t_ind_test_p(predex, predey) < alpha) t_test_WN_rejected = t_test_WN_rejected + 1
    
    ### standardni t-test - ARMA ###
    
    if(cor.test(x, y)$p.value < alpha) standard_t_test_ARMA_rejected = standard_t_test_ARMA_rejected + 1
    
    ### standardni t-test - bily sum ###
    
    if(cor.test(predex, predey)$p.value < alpha) standard_t_test_WN_rejected = standard_t_test_WN_rejected + 1
    
    ### test distancni kovarianci ###
    
    if(dist_ind_test_p(predex, predey) < alpha) dist_test_rejected = dist_test_rejected + 1
    
  }
  
  ### Ulozeni informaci a vysledku do data framu ###
  
  results = data.frame( Number_of_tests,
                        alpha,
                        N,
                        Q,
                        s1,
                        s2,
                        rho,
                        M,
                        mx,
                        nx,
                        paste(ax, collapse = ", "),
                        paste(bx, collapse = ", "),
                        my,
                        ny,
                        paste(ay, collapse = ", "),
                        paste(by, collapse = ", "),
                        haugh_test_v1_rejected,
                        haugh_test_v2_rejected,
                        t_test_ARMA_rejected,
                        t_test_WN_rejected,
                        standard_t_test_ARMA_rejected,
                        standard_t_test_WN_rejected,
                        dist_test_rejected,
                        model_estimation_x_right,
                        model_estimation_y_right)
  return(results)
}

transform_results <- function(results){ # funkce na transformaci data framu z funkce Simulation
  rhos = rep(-10:10/10, 7)
  Names = c(rep('H1', 21),
            rep('H2', 21),
            rep('Mt-ARMA', 21),
            rep('Mt-WN', 21),
            rep('DCT', 21),
            rep('St-ARMA', 21),
            rep('St-WN', 21))
  Names = factor(Names, levels = c('H1','H2','Mt-ARMA','Mt-WN','DCT','St-ARMA','St-WN'))
  No_obs = c(rep(results$N[1], 7*21))
  res = c(results$haugh_test_v1_rejected,
          results$haugh_test_v2_rejected,
          results$t_test_ARMA_rejected,
          results$t_test_WN_rejected,
          results$dist_test_rejected,
          results$standard_t_test_ARMA_rejected,
          results$standard_t_test_WN_rejected)
  transformed = data.frame(rhos, Names, res, No_obs)  
  transformed$res = transformed$res/results$Number_of_tests[1]
  return(transformed)
}

#####################################################
#####  Priklad 1 - velmi silne korelovane rady  #####
#####################################################

set.seed(seed)
N <- 50                                               # pocet pozorovani
Q <- 10                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu
rho <- 0.8                                            # korelace bilych sumu

### Generovani rad - bily sum ###

sigma <- matrix(c(s1^2, s1 * s2 * rho, s1 * s2 * rho, s2^2), 2)   # kovariancni matice
e <- mvrnorm((N+Q), mu, sigma)                        # vygenerovani bilych sumu za pomoci funkce mvrnorm z balicku MASS
colnames(e) <- c("ex", "ey")
ex <- e[ , c("ex")]                                   # bily sum prislusny rade x
ey <- e[ , c("ey")]                                   # bily sum prislusny rade y

### Generovani rad - x - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 2
ax = c(-0.95)
bx = c(-0.1, 0.6)

x = generate_ARMA(N, mx, nx, ax, bx, ex, Q)

### Generovani rad - y - ARMA(my, ny) s parametry ay a by ###

my = 2
ny = 1
ay = c(0.3, -0.2)
by = c(-0.7)

y = generate_ARMA(N, my, ny, ay, by, ey, Q)

### vizualizace rad x a y ###

Index = rep(c(1:N), 2)
Hodnoty = c(x, y)
Proměnná = c(rep("X", N), rep("Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-8, 8) + theme_light()

### vizualizace bilych sumu ex a ey ###

Index = rep(c(1:(N+Q)), 2)
Hodnoty = c(ex, ey)
Proměnná = c(rep("E_X", N+Q), rep("E_Y", N+Q))

data_plot_WN = data.frame(Index, Proměnná, Hodnoty)

ggplot(data_plot_WN, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-3, 3) + theme_light()

########################################
#####  Priklad 2 - nezavisle rady  #####
########################################

set.seed(seed)
N <- 50                                               # pocet pozorovani
Q <- 10                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu
rho <- 0                                              # korelace bilych sumu

### Generovani rad - bily sum ###

sigma <- matrix(c(s1^2, s1 * s2 * rho, s1 * s2 * rho, s2^2), 2)   # kovariancni matice
e <- mvrnorm((N+Q), mu, sigma)                        # vygenerovani bilych sumu za pomoci funkce mvrnorm z balicku MASS
colnames(e) <- c("ex", "ey")
ex <- e[ , c("ex")]                                   # bily sum prislusny rade x
ey <- e[ , c("ey")]                                   # bily sum prislusny rade y

### Generovani rad - x - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 2
ax = c(0.3)
bx = c(-0.7, 0.6)

x = generate_ARMA(N, mx, nx, ax, bx, ex, Q)

### Generovani rad - y - ARMA(my, ny) s parametry ay a by ###

my = 2
ny = 1
ay = c(0.3, -0.2)
by = c(-0.7)

y = generate_ARMA(N, my, ny, ay, by, ey, Q)

### vizualizace rad x a y ###

Index = rep(c(1:N), 2)
Hodnoty = c(x, y)
Proměnná = c(rep("X", N), rep("Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-8, 8) + theme_light()

### vizualizace bilych sumu ex a ey ###

Index = rep(c(1:(N+Q)),2)
Hodnoty = c(ex,ey)
Proměnná = c(rep("E_X",N+Q),rep("E_Y",N+Q))

data_plot_WN = data.frame(Index, Proměnná, Hodnoty)

ggplot(data_plot_WN, aes(x=Index, y = Hodnoty, col= Proměnná))+geom_line(aes(color = Proměnná), lwd = 1)+geom_point(aes(color = Proměnná), size = 2) + ylim(-4, 4) + theme_light()

##########################################
#####  jednoduchy model na simulaci  #####
##########################################

set.seed(seed)
N <- 50                                               # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu
rho <- 0.5                                            # korelace bilych sumu

### Generovani rad - bily sum ###

sigma <- matrix(c(s1^2, s1 * s2 * rho, s1 * s2 * rho, s2^2), 2)   # kovariancni matice
e <- mvrnorm((N+Q), mu, sigma)                        # vygenerovani bilych sumu za pomoci funkce mvrnorm z balicku MASS
colnames(e) <- c("ex", "ey")
ex <- e[ , c("ex")]                                   # bily sum prislusny rade x
ey <- e[ , c("ey")]                                   # bily sum prislusny rade y

### Generovani rad - x - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 0
ax = c(0.8)
bx = 0

x = generate_ARMA(N, mx, nx, ax, bx, ex, Q)

### Generovani rad - y - ARMA(my, ny) s parametry ay a by ###

my = 0
ny = 1
ay = 0
by = c(-0.4)

y = generate_ARMA(N, my, ny, ay, by, ey, Q)

### vizualizace rad x a y ###

Index = rep(c(1:N), 2)
Hodnoty = c(x, y)
Proměnná = c(rep("X", N), rep("Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-5, 5) + theme_light()

### odhad modelu ARMA ###

model_x = find_ARMA_model(x, 3, 3)
summary(model_x)

model_y = find_ARMA_model(y, 3, 3)
summary(model_y)

x_res = model_x$residuals[1:(length(x))]
y_res = model_y$residuals[1:(length(y))]

### vizualizace odhadnutych bilych sumu ###

Index = rep(c(1:N), 2)
Hodnoty = c(x_res, y_res)
Proměnná = c(rep("E_X", N), rep("E_Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-3, 3) + theme_light()

### vyberova krizova korelace - x a y ###

x_y_CC = c()
Lag = c()

for(i in -(length(x) - 1):(length(x) - 1)){    # vypocet vyberove korelace pro vsechny mozne body
  x_y_CC = c(x_y_CC, cross_correlation(x, y, i, WN = FALSE))
  Lag = c(Lag, i)
}
cross_correlation(x, y, 0, WN = FALSE)
x_y_CC_data = data.frame(x_y_CC, Lag)

x_y_CC_graph<- ggplot(x_y_CC_data, aes(x = Lag, y = x_y_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.6, 0.6)
x_y_CC_graph # 5x7 landscape

### realna krizova korelace - x a y ###

R_X_Y = c()
for(k in -(N - 1):0){
  R_X_Y = c(R_X_Y, (-0.8)^abs(k) - 0.4*(-0.8)^(abs(k) + 1))
}
R_X_Y = c(R_X_Y, -0.4)
R_X_Y = c(R_X_Y, rep(0, N - 2))
R_X_Y = R_X_Y*rho  # krizova kovariancni funkce

r_X_Y = R_X_Y/(sqrt(25/9 * 29/25)) # krizova korelacni funkce

x_y_CC_data = data.frame(Value = c(x_y_CC, r_X_Y), Typ = c(rep('Výběrová', 2*N - 1), rep('Reálná', 2*N - 1)), Lag = rep(Lag, 2))

x_y_CC_graph<- ggplot(x_y_CC_data, aes(x = Lag, y = Value, color = Typ)) + geom_line(lwd = 0.5, aes(color = Typ)) + geom_point(aes(color = Typ), size = 1) + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.6, 0.6)
x_y_CC_graph # 5x7 landscape

### vyberova krizova korelace - odhadnute bile sumy ###

x_y_WN_CC = c()
Lag = c()

for(i in -(length(x_res) - 1):(length(x_res) - 1)){    # vypocet vyberove korelace pro vsechny mozne body
  x_y_WN_CC = c(x_y_WN_CC, cross_correlation(x_res, y_res, i, WN = TRUE))
  Lag = c(Lag, i)
}

cross_correlation(x_res, y_res, 0, WN = TRUE)
x_y_WN_CC_data = data.frame(x_y_WN_CC, Lag)

x_y_WN_CC_graph<- ggplot(x_y_WN_CC_data, aes(x = Lag, y = x_y_WN_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.6, 0.6)
x_y_WN_CC_graph # 5x7 landscape


#########################################
#####  komplexni model na simulaci  #####
#########################################

set.seed(seed)
N <- 50                                               # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu
rho <- -0.5                                           # korelace bilych sumu

### Generovani rad - bily sum ###

sigma <- matrix(c(s1^2, s1 * s2 * rho, s1 * s2 * rho, s2^2), 2)   # kovariancni matice
e <- mvrnorm((N+Q), mu, sigma)                        # vygenerovani bilych sumu za pomoci funkce mvrnorm z balicku MASS
colnames(e) <- c("ex", "ey")
ex <- e[ , c("ex")]                                   # bily sum prislusny rade x
ey <- e[ , c("ey")]                                   # bily sum prislusny rade y

### Generovani rad - x - ARMA(mx, nx) s parametry ax a bx ###

mx = 2
nx = 2
ax = c(-0.95, 0.1)
bx = c(-0.1, -0.5)

x = generate_ARMA(N, mx, nx, ax, bx, ex, Q)

### Generovani rad - y - ARMA(my, ny) s parametry ay a by ###

my = 1
ny = 3
ay = c(0.5)
by = c(0.6, 0.2, -0.5)

y = generate_ARMA(N, my, ny, ay, by, ey, Q)

### vizualizace rad x a y ###

Index = rep(c(1:N), 2)
Hodnoty = c(x, y)
Proměnná = c(rep("X", N), rep("Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-5, 5) + theme_light()

### odhad modelu ARMA ###

model_x = find_ARMA_model(x, 3, 3)
summary(model_x)

model_y = find_ARMA_model(y, 3, 3)
summary(model_y)

x_res = model_x$residuals[1:(length(x))]
y_res = model_y$residuals[1:(length(y))]

### vizualizace odhadnutych bilych sumu ###

Index = rep(c(1:N), 2)
Hodnoty = c(x_res, y_res)
Proměnná = c(rep("E_X", N), rep("E_Y", N))

data_plot = data.frame(Index, Hodnoty, Proměnná)

ggplot(data_plot, aes(x=Index, y = Hodnoty, col= Proměnná)) + geom_line(aes(color = Proměnná), lwd = 1) + geom_point(aes(color = Proměnná), size = 2) + ylim(-3, 3) + theme_light()

### vyberova krizova korelace - x a y ###

x_y_CC = c()
Lag = c()

for(i in -(length(x) - 1):(length(x) - 1)){   # vypocet vyberove korelace pro vsechny mozne body
  x_y_CC = c(x_y_CC, cross_correlation(x, y, i, WN = FALSE))
  Lag = c(Lag, i)
}

cross_correlation(x, y, 0, WN = FALSE)
x_y_CC_data = data.frame(x_y_CC, Lag)

x_y_CC_graph<- ggplot(x_y_CC_data, aes(x = Lag, y = x_y_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.5, 0.5)
x_y_CC_graph # 5x7 landscape

### realna krizova korelace - x a y ###

# faktory pro vyjadreni x jako linearniho procesu

g_X = c(1)

for(k in (1:200)){
  g_X = c(g_X, 40^(-k) * ((19 + sqrt(201))^(k) * (-40 + 3 * sqrt(201)) + (19 - sqrt(201))^(k) * (40 + 3 * sqrt(201)))/sqrt(201))
}

# faktory pro vyjadreni y jako linearniho procesu
# zdroj: https://www.wolframalpha.com/input?i=200+terms+of+taylor+series+of+%281%2B0.5x%29%2F%281%2B0.6x%2B0.2x%5E2-0.5x%5E3%29+at+x%3D0
g_Y =   c(1, -0.1, -0.14, 0.604, -0.3844, 0.03984, 0.354976, -0.413154, 0.196817, 0.142029, -0.331157, 0.268697, -0.0239726, -0.204935, 0.262104, -0.128262, -0.077931,
          0.203463, -0.170622, 0.0227153, 0.122227, -0.16319, 0.0848265, 0.0428555, -0.124274, 0.108406, -0.0187613, -0.0725614, 0.101492, -0.0557637, -0.0231209, 0.0757714,
          -0.0687205, 0.0145176, 0.0429193, -0.0630153, 0.0364842, 0.0121722, -0.0461078, 0.0434723, -0.0107757, -0.0252829, 0.0390611, -0.0237679, -0.00619293, 0.0279999,
          -0.0274453, 0.00777074, 0.0148266, -0.0241727, 0.0154237, 0.00299361, -0.0169673, 0.0172935, -0.00548584, -0.00865083, 0.0149344, -0.0099734, -0.00132826, 0.0102588,
          -0.0108764, 0.00380991, 0.00501874, -0.0092114, 0.00642805, 0.000494821, -0.00618821, 0.00682798, -0.00261174, -0.00289266, 0.00567193, -0.0041305, -0.000102415,
          0.00372352, -0.00427888, 0.00177141, 0.00165468, -0.00348653, 0.00264669, -0.0000633651, -0.00223458, 0.00267677, -0.00119083, -0.00093815, 0.00213944, -0.00169145,
          0.000117905, 0.00133727, -0.00167166, 0.000794498, 0.000526267, -0.00131049, 0.00107829, -0.000121743, -0.000797859, 0.00104221, -0.000526625, -0.000291396, 0.000801267,
          -0.000685794, 0.000105525, 0.000474478, -0.000648688, 0.00034708, 0.000158729, -0.000488997, 0.000435192, -0.0000839518, -0.000281166, 0.000403086, -0.000227594,
          -0.0000846436, 0.000297848, -0.000275577, 0.000063455, 0.000165967, -0.00025006, 0.00014857, 0.0000438532, -0.000181056, 0.000174148, -0.0000463509, -0.0000975469,
          0.000154872, -0.0000965894, -0.0000217942, 0.000109831, -0.000109834, 0.0000330373, 0.0000570597, -0.0000957604, 0.0000625629, 0.0000101442, -0.0000664793, 0.0000691402,
          -0.0000231162, -0.000033198, 0.0000591121, -0.0000403858, -4.18996e-06, 0.0000401472, -0.0000434432, 0.0000159415, 0.0000191973, -0.0000364283, 0.0000259883, 1.29137e-06,
          -0.0000241866, 0.0000272478, -0.0000108657, -0.0000110235, 0.0000224111, -0.0000166748, 1.09401e-08, 0.000014534, -0.00001706, 7.33467e-06, 6.27818e-06, -0.0000137638,
          0.00001067, -5.10142e-07, -8.70983e-06, 0.0000106629, -4.91086e-06, -3.54098e-06, 8.43823e-06, -6.81017e-06, 6.27964e-07, 5.20437e-06, -6.6533e-06, 3.26509e-06,
          1.97379e-06, -5.16394e-06, 4.33615e-06, -5.82006e-07, -3.1e-06, 4.14447e-06, -2.15769e-06, -1.08428e-06, 3.15434e-06, -2.75459e-06, 4.79747e-07, 1.84024e-06,
          -2.57739e-06, 1.41826e-06, 5.84643e-07, -1.92313e-06, 1.74608e-06, -3.70701e-07, -1.08836e-06, 1.6002e-06, -9.27797e-07, -3.07543e-07, 1.17018e-06, -1.1045e-06,
          2.74892e-07, 6.41057e-07, -9.91863e-07, 6.04352e-07, 1.5629e-07, -7.10576e-07)


R_X_Y = c()
for(k in -(N - 1):(N - 1)){   # vypocet realne kovariance mezi x a y pro vsechny mozne body
  if(k>=0){
    g_X_C = g_X[1:(length(g_X) - k)]
    g_Y_C = g_Y[(1 + k):(length(g_Y))]
  }else{
    g_X_C = g_X[(-k + 1):(length(g_X))]
    g_Y_C = g_Y[1:((length(g_Y)) + k)]
  }
  R_X_Y = c(R_X_Y, sum(g_X_C * g_Y_C))
}

R_X_Y = R_X_Y*rho  # realna kovariance mezi x a y

r_X_Y = R_X_Y/(sqrt(sum(g_X^2)) * sqrt(sum(g_Y^2)))  # realna korelace mezi x a y

x_y_CC_data = data.frame(Value = c(x_y_CC, r_X_Y), Typ = c(rep('Výběrová', 2*N - 1), rep('Reálná', 2*N - 1)), Lag = rep(Lag, 2))

x_y_CC_graph<- ggplot(x_y_CC_data, aes(x = Lag, y = Value, color = Typ)) + geom_line(lwd = 0.5, aes(color = Typ)) + geom_point(aes(color = Typ), size = 1) + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.5, 0.5)
x_y_CC_graph # 5x7 landscape


### vyberova krizova korelace - odhadnute bile sumy ###

x_y_WN_CC = c()
Lag = c()

for(i in -(length(x_res) - 1):(length(x_res) - 1)){
  x_y_WN_CC = c(x_y_WN_CC, cross_correlation(x_res, y_res, i, WN = TRUE))
  Lag = c(Lag, i)
}

cross_correlation(x_res, y_res, 0, WN = TRUE)
x_y_WN_CC_data = data.frame(x_y_WN_CC, Lag)

x_y_WN_CC_graph<- ggplot(x_y_WN_CC_data, aes(x = Lag, y = x_y_WN_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.5, 0.5)
x_y_WN_CC_graph # 5x7 landscape

############################################################
#####  Nastaveni simulace 1 - kratky jednoduchy model  #####
############################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 50                                               # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu
                                               
M = 5                                                 # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 0
ax = c(0.8)
bx = 0

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 0
ny = 1
ay = 0
by = c(-0.4)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results1 = results

####################################################################
#####  Nastaveni simulace 2 - stredne dlouhy jednoduchy model  #####
####################################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 100                                              # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu

M = 10                                                # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 0
ax = c(0.8)
bx = 0

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 0
ny = 1
ay = 0
by = c(-0.4)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results2 = results

############################################################
#####  Nastaveni simulace 3 - dlouhy jednoduchy model  #####
############################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 200                                              # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu

M = 20                                                # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 1
nx = 0
ax = c(0.8)
bx = 0

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 0
ny = 1
ay = 0
by = c(-0.4)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results3 = results

###########################################################
#####  Nastaveni simulace 4 - kratky komplexni model  #####
###########################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 50                                               # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu

M = 5                                                # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 2
nx = 2
ax = c(-0.95, 0.1)
bx = c(-0.1, -0.5)

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 1
ny = 3
ay = c(0.5)
by = c(0.6, 0.2, -0.5)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results4 = results

###################################################################
#####  Nastaveni simulace 5 - stredne dlouhy komplexni model  #####
###################################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 100                                              # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu

M = 10                                                # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 2
nx = 2
ax = c(-0.95, 0.1)
bx = c(-0.1, -0.5)

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 1
ny = 3
ay = c(0.5)
by = c(0.6, 0.2, -0.5)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results5 = results

###########################################################
#####  Nastaveni simulace 6 - dlouhy komplexni model  #####
###########################################################

set.seed(seed)
Number_of_tests = 10000                               # pocet simulaci
alpha = 0.05

N <- 200                                              # pocet pozorovani
Q <- 30                                               # pocet pozorovani vygenerovanych navic, ktere se zahodi
mu <- c(0, 0)                                         # vektor strednich hodnot bilych sumu
s1 <- 1                                               # smerodatna odchylka prvniho bileho sumu
s2 <- 1                                               # smerodatna odchylka druheho bileho sumu

M = 20                                                # Haugh test M

### nastaveni rady X - ARMA(mx, nx) s parametry ax a bx ###

mx = 2
nx = 2
ax = c(-0.95, 0.1)
bx = c(-0.1, -0.5)

### nastaveni rady Y - ARMA(my, ny) s parametry ay a by ###

my = 1
ny = 3
ay = c(0.5)
by = c(0.6, 0.2, -0.5)

### simulace ###

step = -10:10/10

for(i in step){  # korelace i
  t = Sys.time()
  print(t)
  print(i)
  if(i == -1){
    results = Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by)
  }else{ 
    results = rbind(results, Simulation(Number_of_tests, alpha, N, Q, mu, s1, s2, i, M, mx, nx, ax, bx, my, ny, ay, by))
  }
}

results6 = results

##########################################
#####  Graficke zpracovani vysledku  #####
##########################################

transformed_1 = transform_results(results1)
transformed_2 = transform_results(results2)
transformed_3 = transform_results(results3)

combined_results_1 = rbind(transformed_1, transformed_2, transformed_3)
combined_results_1$No_obs = factor(combined_results_1$No_obs, levels = c("50", "100", "200"))

combined_results_graph_1<-ggplot(combined_results_1, aes(x = rhos, y = res, colour = No_obs)) + geom_point(size = 4) + facet_wrap(~Names, ncol = 2) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Počet pozorování") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

combined_results_graph_1 # 8x12 portrait

transformed_4 = transform_results(results4)
transformed_5 = transform_results(results5)
transformed_6 = transform_results(results6)

combined_results_2 = rbind(transformed_4, transformed_5, transformed_6)
combined_results_2$No_obs = factor(combined_results_2$No_obs, levels = c("50", "100", "200"))

combined_results_graph_2<-ggplot(combined_results_2, aes(x = rhos, y = res, colour = No_obs)) + geom_point(size = 4) + facet_wrap(~Names, ncol = 2) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Počet pozorování") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

combined_results_graph_2 # 8x12 portrait


results_graph_1<-ggplot(transformed_1, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_1 # 6x8 landscape
 
results_graph_2<-ggplot(transformed_2, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_2 # 6x8 landscape

results_graph_3<-ggplot(transformed_3, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_3 # 6x8 landscape

results_graph_4<-ggplot(transformed_4, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_4 # 6x8 landscape

results_graph_5<-ggplot(transformed_5, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_5 # 6x8 landscape

results_graph_6<-ggplot(transformed_6, aes(x = rhos, y = res, colour = Names))+ geom_point(size = 4) + labs(x = "Korelace mezi bílými šumy", y = "Podíl zamítnutých testů", color = "Název testu") + geom_hline(yintercept = 0.05, linetype = "dashed", color = "red")

results_graph_6 # 6x8 landscape

########################################
#####  Priklad na realnych datech  #####
########################################

### nacteni a uprava dat ###

# zdroj pro S&P 500: https://www.wsj.com/market-data/quotes/index/SPX/historical-prices
# zdroj pro FTSE 100: https://www.wsj.com/market-data/quotes/index/UK/FTSE%20UK/UKX/historical-prices
# Data v obdobi 1. 6. 2021 - 31. 5. 2023

SandP500 <- read.csv(".\\S&P500.csv", header=TRUE, stringsAsFactors=FALSE)
FTSE100  <- read.csv(".\\ftse100.csv", header=TRUE, stringsAsFactors=FALSE)

intersection = inner_join(SandP500, FTSE100, by="Date")   # vybrani dnu, kdy se obchodovalo s obema indexy

SandP500_edit = data.frame('Date' = as.Date(intersection$Date, format = "%m/%d/%y"), 'Value' = intersection$Close.x)
FTSE100_edit = data.frame('Date' = as.Date(intersection$Date, format = "%m/%d/%y"),  'Value' = intersection$Close.y)
SandP500_edit = SandP500_edit[order(nrow(SandP500_edit):1),]
FTSE100_edit = FTSE100_edit[order(nrow(FTSE100_edit):1),]

### vyberove statistiky ###

print(summary(SandP500_edit$Value), digits = 6 )
print(var(SandP500_edit$Value), digits = 7 )
print(summary(FTSE100_edit$Value), digits = 6 )
print(var(FTSE100_edit$Value), digits = 7 )

### vizualizace cen indexu ###

SandP500_graph <- ggplot(SandP500_edit, aes(x = Date, y = Value)) + geom_line(lwd = 0.5, color = "#F8766D") + theme_light() + labs(x = "Datum", y = "Hodnota") + scale_x_date(date_labels = "%m-%y",date_breaks  ="3 month",date_minor_breaks = "1 month") + theme(axis.text.x = element_text(angle = 45, hjust = 1))
SandP500_graph # 6x8 landscape

FTSE100_graph <- ggplot(FTSE100_edit, aes(x = Date, y = Value)) + geom_line(lwd = 0.5, color = "#00BFC4") + theme_light() + labs(x = "Datum", y = "Hodnota") + scale_x_date(date_labels = "%m-%y",date_breaks  ="3 month",date_minor_breaks = "1 month") + theme(axis.text.x = element_text(angle = 45, hjust = 1))
FTSE100_graph # 6x8 landscape

### odhady modelu ARMA ###

modelSP500 = find_ARMA_model_with_mean(SandP500_edit$Value, 3, 3)
summary(modelSP500)

modelFTSE100 = find_ARMA_model_with_mean(FTSE100_edit$Value, 3, 3)
summary(modelFTSE100)

SP500Res = modelSP500$residuals[1:(length(SandP500_edit$Value))]
FTSE100Res = modelFTSE100$residuals[1:(length(FTSE100_edit$Value))]

summary(SP500Res)
summary(FTSE100Res)

### vizualizace odhadhu bilych sumu ###

SP500Resdata = data.frame(Date = SandP500_edit$Date, Value = SP500Res)
FTSE100Resdata = data.frame(Date = FTSE100_edit$Date, Value = FTSE100Res)

SandP500Res_graph <- ggplot(SP500Resdata, aes(x = Date, y = Value)) + geom_line(lwd = 0.5, color = "#F8766D") + theme_light() + labs(x = "Datum", y = "Hodnota") + scale_x_date(date_labels = "%m-%y",date_breaks  ="3 month",date_minor_breaks = "1 month") + theme(axis.text.x = element_text(angle = 45, hjust = 1))
SandP500Res_graph # 6x8 landscape

FTSE100Res_graph <- ggplot(FTSE100Resdata, aes(x = Date, y = Value)) + geom_line(lwd = 0.5, color = "#00BFC4") + theme_light() + labs(x = "Datum", y = "Hodnota") + scale_x_date(date_labels = "%m-%y",date_breaks  ="3 month",date_minor_breaks = "1 month") + theme(axis.text.x = element_text(angle = 45, hjust = 1))
FTSE100Res_graph # 6x8 landscape

combined_res = data.frame(Date = rep(SandP500_edit$Date, 2), Index = c(rep('S&P 500', length(SP500Res)),rep('FTSE 100', length(FTSE100Res))), Value = c(SP500Res, FTSE100Res))

combined_res_graph <- ggplot(combined_res, aes(x = Date, y = Value, col = Index)) + geom_line(lwd = 0.5, aes(color = Index)) + theme_light() + labs(x = "Datum", y = "Hodnota") + scale_x_date(date_labels = "%m-%y",date_breaks  ="3 month",date_minor_breaks = "1 month") + theme(axis.text.x = element_text(angle = 45, hjust = 1))
combined_res_graph

### testovani nezavislosti ###

M = 40

### Haugh test - v1 ###

Haugh_ind_test_p(SP500Res, FTSE100Res, M, better_approx = FALSE)

### Haugh test - v2 ###

Haugh_ind_test_p(SP500Res, FTSE100Res, M, better_approx = TRUE)

### modifikovany t-test - ARMA ###

compute_N_hat(SandP500_edit$Value, FTSE100_edit$Value)
t_ind_test_p(SandP500_edit$Value, FTSE100_edit$Value)

### modifikovany t-test - odhady bilych sumu ###

compute_N_hat(SP500Res, FTSE100Res)
t_ind_test_p(SP500Res, FTSE100Res)

### standardni t-test - ARMA ###

cor(SandP500_edit$Value, FTSE100_edit$Value)
cor.test(SandP500_edit$Value, FTSE100_edit$Value)$p.value

### standardni t-test - odhady bilch sumu ###

cor(SP500Res, FTSE100Res)
cor.test(SP500Res, FTSE100Res)$p.value

### test distancni kovarianci ###

dist_corr(SP500Res, FTSE100Res)
dist_ind_test_p(SP500Res, FTSE100Res)


### vyberova krizova korelace mezi cenami indexu ###

SP_FTSE_CC = c()
Lag = c()

for(i in -(length(SandP500_edit$Value) - 1):(length(SandP500_edit$Value) - 1)){
  SP_FTSE_CC = c(SP_FTSE_CC, cross_correlation(SandP500_edit$Value, FTSE100_edit$Value, i, WN= FALSE))
  Lag = c(Lag, i)
}

SP_FTSE_CC_data = data.frame(SP_FTSE_CC, Lag)

SP_FTSE_CC_graph<- ggplot(SP_FTSE_CC_data, aes(x = Lag, y = SP_FTSE_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.4, 0.4)
SP_FTSE_CC_graph # 5x7 landscape

### vyberova krizova korelace mezi odhady bilych sumu ###

SP_FTSE_WN_CC = c()
Lag = c()

for(i in -(length(SP500Res) - 1):(length(SP500Res) - 1)){
  SP_FTSE_WN_CC = c(SP_FTSE_WN_CC, cross_correlation(SP500Res, FTSE100Res, i, WN = TRUE))
  Lag = c(Lag, i)
}

SP_FTSE_WN_CC_data = data.frame(SP_FTSE_WN_CC, Lag)

SP_FTSE_WN_CC_graph<- ggplot(SP_FTSE_WN_CC_data, aes(x = Lag, y = SP_FTSE_WN_CC)) + geom_line(lwd = 0.5, color = "#000000") + theme_light() + labs(x = "k", y = "Hodnota křížové korelace") + ylim(-0.5, 0.5)
SP_FTSE_WN_CC_graph # 5x7 landscape

### Overeni predpokladu ###

# normalita
shapiro.test(SP500Res)
shapiro.test(FTSE100Res)

# nekorelovanost
Box.test(SP500Res, lag = 50, type = c("Ljung-Box"))
Box.test(FTSE100Res, lag = 50, type = c("Ljung-Box"))