Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(expm)
- ###### polynomial ######
- p <- 4
- mt <- rep(1, p)
- G <- diag(p); G[1:(p-1), 2:p] <- G[1:(p-1), 2:p] + diag(p-1);
- F <- c(1, rep(0, p-1))
- print(G); print(mt); print(F)
- k <- 5
- G %^% k
- F * ((G %^% k) %*% mt)
- ############ seasonal ##########
- mt <- rep(1, 2); F <- c(1,0); w <- pi/5; r <- 0.9 # set r to 1 and it won't damp
- G <- r*matrix(c(cos(w), sin(w), -sin(w), cos(w)), nrow =2, byrow = TRUE)
- k <- 50
- y <- sapply(1:k, function(i) sum(F * ((G %^% i) %*% mt)))
- plot(y, type = 'l')
- ######## combo seosonal #########
- mt1 <- rep(1, 2); F1 <- c(1,0); w1 <- pi/5; r1 <- 0.95
- G1 <- r1*matrix(c(cos(w1), sin(w1), -sin(w1), cos(w1)), nrow =2, byrow = TRUE)
- y1 <- sapply(1:k, function(i) sum(F1 * ((G1 %^% i) %*% mt1)))
- mt2 <- rep(2, 2); F2 <- c(1,0); w2 <- pi/10; r2 <- 0.95
- G2 <- r2*matrix(c(cos(w2), sin(w2), -sin(w2), cos(w2)), nrow =2, byrow = TRUE)
- y2 <- sapply(1:k, function(i) sum(F2 * ((G2 %^% i) %*% mt2)))
- mt3 <- rep(2, 2); F3 <- c(1,0); w3 <- pi/20; r3 <- 0.95
- G3 <- r3*matrix(c(cos(w3), sin(w3), -sin(w3), cos(w3)), nrow =2, byrow = TRUE)
- y3 <- sapply(1:k, function(i) sum(F3 * ((G3 %^% i) %*% mt3)))
- mt <- c(mt1, mt2, mt3); F <- c(F1, F2, F3); G <- as.matrix(bdiag(G1, G2, G3))
- y <- sapply(1:k, function(i) sum(F * ((G %^% i) %*% mt)))
- par(mfrow=c(2,1), mar = c(2, 2, 2, 2))
- plot(y1, ylim = c(-4,4), type = 'l', col = "blue"); lines(y2, col = "red"); lines(y3, col = "green")
- plot(y, type = 'l')
Add Comment
Please, Sign In to add comment