Guest User

Untitled

a guest
Oct 24th, 2018
126
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.39 KB | None | 0 0
  1. library(expm)
  2.  
  3. ###### polynomial ######
  4.  
  5. p <- 4
  6. mt <- rep(1, p)
  7. G <- diag(p); G[1:(p-1), 2:p] <- G[1:(p-1), 2:p] + diag(p-1);
  8. F <- c(1, rep(0, p-1))
  9.  
  10. print(G); print(mt); print(F)
  11.  
  12. k <- 5
  13. G %^% k
  14.  
  15. F * ((G %^% k) %*% mt)
  16.  
  17.  
  18. ############ seasonal ##########
  19.  
  20. mt <- rep(1, 2); F <- c(1,0); w <- pi/5; r <- 0.9 # set r to 1 and it won't damp
  21. G <- r*matrix(c(cos(w), sin(w), -sin(w), cos(w)), nrow =2, byrow = TRUE)
  22.  
  23. k <- 50
  24. y <- sapply(1:k, function(i) sum(F * ((G %^% i) %*% mt)))
  25. plot(y, type = 'l')
  26.  
  27. ######## combo seosonal #########
  28. mt1 <- rep(1, 2); F1 <- c(1,0); w1 <- pi/5; r1 <- 0.95
  29. G1 <- r1*matrix(c(cos(w1), sin(w1), -sin(w1), cos(w1)), nrow =2, byrow = TRUE)
  30. y1 <- sapply(1:k, function(i) sum(F1 * ((G1 %^% i) %*% mt1)))
  31.  
  32. mt2 <- rep(2, 2); F2 <- c(1,0); w2 <- pi/10; r2 <- 0.95
  33. G2 <- r2*matrix(c(cos(w2), sin(w2), -sin(w2), cos(w2)), nrow =2, byrow = TRUE)
  34. y2 <- sapply(1:k, function(i) sum(F2 * ((G2 %^% i) %*% mt2)))
  35.  
  36. mt3 <- rep(2, 2); F3 <- c(1,0); w3 <- pi/20; r3 <- 0.95
  37. G3 <- r3*matrix(c(cos(w3), sin(w3), -sin(w3), cos(w3)), nrow =2, byrow = TRUE)
  38. y3 <- sapply(1:k, function(i) sum(F3 * ((G3 %^% i) %*% mt3)))
  39.  
  40. mt <- c(mt1, mt2, mt3); F <- c(F1, F2, F3); G <- as.matrix(bdiag(G1, G2, G3))
  41. y <- sapply(1:k, function(i) sum(F * ((G %^% i) %*% mt)))
  42.  
  43. par(mfrow=c(2,1), mar = c(2, 2, 2, 2))
  44. plot(y1, ylim = c(-4,4), type = 'l', col = "blue"); lines(y2, col = "red"); lines(y3, col = "green")
  45. plot(y, type = 'l')
Add Comment
Please, Sign In to add comment