install.packages("tsDyn")
library(vars)
library(urca)
library(stargazer)
library(tsDyn)

data <- Dados_Problema12[33:135,2:17]

# padronização
padr <- function(x){
  x = (x - mean(x))/sd(x)
}
data <- as.data.frame(lapply(data, padr))

# 1a diferença // ???
#data[[1]] <-c(NA,diff(data[[1]]))
#data[[3]] <-c(NA,diff(data[[3]]))
#data <- data[-1,]

# VAR
select <- VARselect(data)
reg_var <- VAR(data, p = select$selection[[4]], type = "const")
#stargazer(reg_var$varresult, type = "text")

# PFA
omega_hat <- var(as.matrix(data))
c_hat <- (unlist(lapply(data,mean)))

autovetores <- eigen(omega_hat)$vector
colnames(auto_vetores) <- colnames(data)

autovalores <- as.data.frame(t(eigen(omega_hat)$value)) # auto valores decrescentes
colnames(auto_valores) <- colnames(data)
autovalores <- sort(auto_valores, decreasing = T)

# Variância Explicada // Não tenho certeza
i <- 1:16
var_explicada <- autovalores[i]/sum(autovalores[i])
stargazer(var_explicada, type = 'text',summary.stat="mean") # k = 2

z1 <- as.matrix(data) %*% autovetores[,1]
z2 <- as.matrix(data) %*% autovetores[,2]

# PCA
W <- autovetores[,1:2]
Z <- cbind(z1,z2)
lambda <- t(Z) %*% Z
lambda

pi <- W %*% lambda**(1/2)
D <- diag(omega_hat - pi %*% t(pi))
stargazer(D, type = "text")
View(D)

# Funções para PFA
auto_vec_2 <- function(s){
  eigen(omega_fun(s))$vector[,1:2]
}

lambda_fun <- function(s){
  Z = as.matrix(data) %*% auto_vec_2(s)
  t(Z) %*% Z
}
lambda_fun(1)

pi_fun <- function(s){
  auto_vec_2(s-1) %*% lambda_fun(s-1)**(1/2)
}
pi
omega_fun <- function(s){
  if (s == 0) omega_hat
  else        omega_fun(s-1) - D_fun(s-1)
}

D_fun <- function(s){
  if (s == 0) 0
  else       diag(omega_fun(s) - pi_fun(s-1) %*% t(pi_fun(s-1)))
}

pfa_fun <- function(delta){
  s <- 1
  while(prod(D_fun(s)-D_fun(s-1)) > delta){
    print(s)
    s <- s+1
  }
  invisible(list(pi_fun(s),D_fun(s)))
}