library(tidyverse)

phis <- c(0.6, -0.4, 0.2, 0.2, 0.1)
p <- length(phis)
X <- rnorm(p)

G <- rbind(phis,
           cbind(diag(p-1),
                 rep(0, p-1)))

F <- c(1, rep(0, p-1))

E <- eigen(G)$vectors
A <- eigen(G)$values

# confirm eigendecomp work
round(E %*% diag(A) %*% solve(E), 4) == G

Z <- solve(E) %*% X
F0 <- t(E) %*% F

k <- 50
Y <- data.frame(matrix(NA, k, p))

for (i in 1:p){

  if (Im(A[i]) == 0) {
    Y[, i] <- Re(sapply(1:k, function(t) A[i]^t*F0[i]*Z[i]))
  }

  if (Im(A[i]) != 0) {
    match <- which(Im(A) == -Im(A[i]))
    Y[, i] <- Re(sapply(1:k, function(t) A[i]^t*F0[i]*Z[i])) +
      Re(sapply(1:k, function(t) A[match]^t*F0[match]*Z[match]))
  }

  names(Y)[i] <- paste0("Eig", "_Mod", round(Mod(A[i]), 2), "_Period", round(1/Arg(A[i]), 2))

}

# remove duplicate columns
Y.df <- data.frame((t(unique(t(Y)))))


Y.df <- Y.df %>% mutate(
  t = 1:k,
  Total = Re(sapply(1:k, function(i) t(F0) %*% diag(A)^i %*% Z))
)

Y.df %>%
  gather(Y, val, -t) %>%
  ggplot(aes(t, val, col = Y)) + geom_smooth(se = FALSE, span = 0.05) +
  facet_wrap(~Y, scales = "free", nrow = ncol(Y.df) - 1) +
  theme(legend.position="none") +
  labs(y = "")