install.packages("prophet")
install.packages("prophet")
library(prophet)
install.packages("prophet")
install.packages("prophet")
install.packages("bfast")
library(bfast)
?bfast
t <- ts(start=1,end=4*365,frequency=)
t <- ts(start=1,end=4*365,frequency=1)
t
t <- ts(start=1,end=4*365,frequency=365)
t]
t
t <- ts(1:(4*365),frequency=365)
t
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
beta1 <- c(0,0.1,0.1,0.1)
beta2 <- c(0.5,0,0,0)
X1 <- X[1:(2*365),]
X2 <- X[(2*365+1):(4:365),]
X2 <- X[(2*365+1):(4*365),]
Y1 <- X1*beta1
Y2 <- X2*beta2
Y <- rbind(Y1,Y2)
Y
Y1
Y1 <- X1%*%beta1
Y2 <- X2%*%beta2
Y1
Y <- as.ts(c(Y1,Y2),frequency=365)
Y
Y1 <- X1%*%beta1+rnorm(2*365,0,01)
Y1 <- X1%*%beta1+rnorm(2*365,0,0.1)
Y2 <- X2%*%beta2+rnorm(2*365,0,0.1)
Y <- as.ts(c(Y1,Y2),frequency=365)
plot(Y)
bfast(Y)
length(Y)
head(Y)
Y <- as.ts(c(Y1,Y2),frequency=365)
head(Y)
Y <- ts(c(Y1,Y2),frequency=365)
head(Y)
length(Y)
?ts
ts(Y,start=c(1,1),end=c(4,365))
ts(Y,start=c(1,1),end=c(4,365),frequency=1)
ts(Y,start=c(1,1),end=c(4,365),frequency=365)
Y <- ts(c(Y1,Y2),start=c(1,1),end=c(4,365),frequency=365)
head(Y)
Y
head(Y)
Y <- as.ts(c(Y1,Y2),start=c(1,1),end=c(4,365),frequency=365)
head(Y)
Y <- ts(c(Y1,Y2),start=c(1,1),frequency=365)
head(Y)
Y <- ts(c(Y1,Y2),start=c(1,1),frequency=1)
head(Y)
length(c(Y1,Y2))
length(c(Y1,Y2))/365
Yx <- c(Y1,Y2)
plot(Yx)
ts(Yx)
plot(ts(Yx))
plot(ts(Yx,frequency=365))
head(ts(Yx,frequency=365))
Y <- ts(Yx,deltat=1)
head(Y)
Y <- ts(Yx,frequency=365)
bfast(Y)
mod <- bfast(Y)
mod$jump
abline(v=mod$jump$x)
Y.ts <- Y
mod$output
mod$nobp
mod$Mags
mod$Magnitude
mod$Time
?split
365*2
# extract break times
tchange <- mod$Time
#Segment data
split(Y.ts, cut(Y.ts, c(-Inf, tchange, Inf)))
#Segment data
split(c(Y.ts), cut(Y.ts, c(-Inf, tchange, Inf)))
#Segment data
as.matrix(split(c(Y.ts), cut(Y.ts, c(-Inf, tchange, Inf))))
t <- 1:(365*4)
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
X
X1 <- X[1:1*365,]
X2 <- X[(365+1):(365*2),]
X3 <- X[(365*2+1):(365*4),]
beta1 <- c(0,0.1,0.1,0.1)
beta2 <- c(0.5,0,0,0)
beta3 <- c(0,0.05,0,0)
Y1 <- X1%*%beta1 + rnorm(365,0,0.2)
X1
X1 <- X[1:365,]
X1
Y1 <- X1%*%beta1 + rnorm(365,0,0.2)
Y2 <- X2%*%beta2 + rnorm(365,0,0.2)
Y3 <- X3%*%beta3 + rnorm(365,0,0.2)
Y <- c(Y1,Y2,Y3)
#### Generate some data
t <- 1:(365*4)
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
X1 <- X[1:1*365,]
X2 <- X[(365+1):(365*2),]
X3 <- X[(365*2+1):(365*4),]
beta1 <- c(0,0.1,0.1,0.1)
beta2 <- c(0.5,0,0,0)
beta3 <- c(0,0.05,0,0)
Y1 <- X1%*%beta1 + rnorm(365,0,0.2)
Y2 <- X2%*%beta2 + rnorm(365,0,0.2)
Y3 <- X3%*%beta3 + rnorm(365,0,0.2)
#### Generate some data
t <- 1:(365*4)
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
X1 <- X[1:365,]
X2 <- X[(365+1):(365*2),]
X3 <- X[(365*2+1):(365*4),]
beta1 <- c(0,0.1,0.1,0.1)
beta2 <- c(0.5,0,0,0)
beta3 <- c(0,0.05,0,0)
Y1 <- X1%*%beta1 + rnorm(365,0,0.2)
Y2 <- X2%*%beta2 + rnorm(365,0,0.2)
Y3 <- X3%*%beta3 + rnorm(365,0,0.2)
Y <- c(Y1,Y2,Y3)
#### Generate some data
t <- 1:(365*4)
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
X1 <- X[1:365,]
X2 <- X[(365+1):(365*2),]
X3 <- X[(365*2+1):(365*4),]
beta1 <- matrix(c(0,0.1,0.1,0.1),4,2)
beta2 <- matrix(c(0.5,0,0,0),4,2)
beta3 <- matrix(c(0,0.05,0,0),4,2)
#### Generate some data
t <- 1:(365*4)
X <- cbind(1,sin(2*pi*t/365),cos(2*pi*t/365),t/1000)
X1 <- X[1:365,]
X2 <- X[(365+1):(365*2),]
X3 <- X[(365*2+1):(365*4),]
beta1 <- matrix(c(0,0.1,0.1,0.1),4,2)
beta2 <- matrix(c(0.5,0,0,0),4,2)
beta3 <- matrix(c(0,0.05,0,0),4,2)
Y1 <- X1%*%beta1 + matrix(rnorm(2*365,0,0.2),365,2)
Y2 <- X2%*%beta2 + matrix(rnorm(2*365,0,0.2),365,2)
Y3 <- X3%*%beta3 + matrix(rnorm(2*365,0,0.2),365,2)
Y1 <- X1%*%beta1 + matrix(rnorm(2*365,0,0.2),365,2)
Y2 <- X2%*%beta2 + matrix(rnorm(2*365,0,0.2),365,2)
Y3 <- X3%*%beta3 + matrix(rnorm(4*365,0,0.2),2*365,2)
Y <- rbind(Y1,Y2,Y3)
plot(Y[,1])
plot(Y[,2])
Y1 <- X1%*%beta1 + matrix(rnorm(2*365,0,0.1),365,2)
Y2 <- X2%*%beta2 + matrix(rnorm(2*365,0,0.1),365,2)
Y3 <- X3%*%beta3 + matrix(rnorm(4*365,0,0.1),2*365,2)
Y <- rbind(Y1,Y2,Y3)
plot(Y[,1])
plot(Y[,2])
# get data in form of ts
Y.ts <- ts(c(t(Y)),frequency=2*365)
Y.ts
plot(Y.ts)
beta1 <- matrix(c(0,0.1,0.1,0.1,0,0.1,0.1,0.1),4,2)
beta2 <- matrix(c(0.5,0,0,0,-0.5,0,0,0),4,2)
beta3 <- matrix(c(0,0.05,0,0,0,-0.05,0,0),4,2)
Y1 <- X1%*%beta1 + matrix(rnorm(2*365,0,0.1),365,2)
Y2 <- X2%*%beta2 + matrix(rnorm(2*365,0,0.1),365,2)
Y3 <- X3%*%beta3 + matrix(rnorm(4*365,0,0.1),2*365,2)
Y <- rbind(Y1,Y2,Y3)
plot(Y[,1])
plot(Y[,2])
k <- 4
d <- 2
?bfast
install.packages("bfast")
#library(rjags)
#library(splines2)
library(bfast)
?bfast
# get data in form of ts
Y.ts <- ts(c(t(Y)),frequency=2*365)
Y_vec <- c(t(Y))
X_vec <- do.call(rbind, lapply(
1:nrow(X),
function(x) {
t(cbind(
c(X[x, ], rep(0, k)),
c(rep(0, k), X[x, ])
))
}
))
X_vec
dim(X_vec)
bfast(Y.ts,season="none",xreg=X_vec)
mod <- bfast(Y.ts,season="none",xreg=X_vec)
mod$jump
mod$output
mod$nobp
mod$Time
plot(Y[,1])
abline(549/2)
abline(v=549/2)
mod$output[[1]]
cbind(Y.ts,X_vec)
df <- cbind(Y.ts,X_vec)
head(df)
bpp <- bfastpp(df,order=0)
breakpoints(response ~ xreg,data=bpp)
mod <- breakpoints(response ~ xreg,data=bpp)
bpp
head(bpp)
dim(bpp)
?bfastpp
?breakpoints
mod <- breakpoints(response ~ xreg,data=bpp,breaks=6)
bpp <- bfastpp(df,order=0)
mod <- breakpoints(response ~ xreg,data=bpp,breaks=6)
bpp <- bfastpp(df,order=0)
mod <- breakpoints(response ~ xreg,data=bpp,breaks=3)
mod <- breakpoints(response ~ xreg,data=bpp,breaks=1)
mod <- breakpoints(response ~ xreg-1,data=bpp,breaks=1)
mod$X
dim(mod$X)
mod$breakpoints
plot(Y[,1])
plot(1460/2)
plot(Y[,1])
abline(v=1460/2)
mod <- breakpoints(response ~ xreg-1,data=bpp,breaks=2)
mod$breakpoints
mod <- breakpoints(response ~ xreg-1,data=bpp,breaks=3)
mod$breakpoints
mod <- breakpoints(response ~ xreg-1,data=bpp,breaks="LWZ")
mod$breakpoints
abline(v=mod$breakpoints)
abline(v=mod$breakpoints/2)
plot(Y[,1])
abline(v=mod$breakpoints/2)
mod$breakpoints/2
#Segment data
as.matrix(split(c(Y.ts), cut(Y.ts, c(-Inf, tchange, Inf))))
# extract break times
tchange <- mod$breakpoints/d
#Segment data
as.matrix(split(c(Y.ts), cut(Y.ts, c(-Inf, tchange, Inf))))
tchange
#Segment data
as.matrix(split(c(Y), cut(Y, c(-Inf, tchange, Inf))))
tchange
#Segment data
thing <- as.matrix(split(c(Y), cut(Y, c(-Inf, tchange, Inf))))
thing
head(thing)
thing[1,]
?split
length(thing)
dim(thing)
thing[[1]]
#Segment data
Y.segs <- as.matrix(split(c(Y), cut(Y, c(-Inf, tchange, Inf))))
X.segs <- split(c(X), cut(X, c(-Inf, tchange, Inf)))
X.segs
dim(X)
cut(X, c(-Inf, tchange, Inf))
splie
X.segs <- split((X), cut(X, c(-Inf, tchange, Inf)))
X.segs
cut(X, c(-Inf, tchange, Inf))
X.segs <- split(data.frame(X), cut(X, c(-Inf, tchange, Inf)))
X.segs
dim(data.frame(X))
length(cut(X, c(-Inf, tchange, Inf)))
X.segs <- split(data.frame(X), cut(Y, c(-Inf, tchange, Inf)))
length(cut(Y, c(-Inf, tchange, Inf)))
dim(Y)
#Segment data
Y.segs <- as.matrix(split(c(Y), cut(Y[,1], c(-Inf, tchange, Inf))))
X.segs <- split(data.frame(X), cut(Y[,1], c(-Inf, tchange, Inf)))
Y.segs
#Segment data
Y.segs <- as.matrix(split((Y), cut(Y[,1], c(-Inf, tchange, Inf))))
Y.segs
#Segment data
Y.segs <- (split((Y), cut(Y[,1], c(-Inf, tchange, Inf))))
Y.segs
#Segment data
Y.segs <- (split(data.frame(Y), cut(Y[,1], c(-Inf, tchange, Inf))))
Y.segs
cut(Y[,1], c(-Inf, tchange, Inf)))
cut(Y[,1], c(-Inf, tchange, Inf))
length(cut(Y[,1], c(-Inf, tchange, Inf))))
length(cut(Y[,1], c(-Inf, tchange, Inf)))
data.frame(Y)
#Segment data
Y.segs <- (split(data.frame(Y), cut(Y[,1], c(-Inf, tchange, Inf))))
Y.segs
dim(Y.segs[[1]])
?cut
#Segment data
Y.segs <- (split(data.frame(Y), cut(1:nrow(Y), c(-Inf, tchange, Inf))))
Y.segs
length(Y.segs)
head(Y.segs[[1]])
head(Y.segs[[2]])
#Segment data
Y.segs <- (split(data.frame(Y), cut(1:nrow(Y), c(-Inf, tchange, Inf))))
X.segs <- split(data.frame(X), cut(1:nrow(Y), c(-Inf, tchange, Inf)))
Y.seg <- Y.segs[[1]]
X.seg <- X.segs[[2]]
nrow(Y.seg) > (365)
nrow(Y.seg)
nrow(Y.seg) >= (365)
Bhat <- matrix(0, k, d)
resids <- matrix(0, n, d)
for (i in 1:d) {
fm1 <- lm(Y.seg[, i] ~ -1 + X.seg)
# allresids <- fm1$residuals
# allresids[abs(allresids)>1] <- NA
resids[, i] <- fm1$residuals
Bhat[, i] <- fm1$coefficients
}
dim(X.seg)
class(list)
class(X.seg)
fm1 <- lm(Y.seg[, i] ~ -1 + as.matrix(X.seg))
fm1
# allresids <- fm1$residuals
# allresids[abs(allresids)>1] <- NA
resids[, i] <- fm1$residuals
Bhat[, i] <- fm1$coefficients
Sigmahat <- cov(resids)
par_vals <- list(Sigmahat = Sigmahat, Bhat = Bhat)
Bhat <- matrix(0, k, d)
resids <- matrix(0, n, d)
n <- nrow(Y)
Bhat <- matrix(0, k, d)
resids <- matrix(0, n, d)
for (i in 1:d) {
fm1 <- lm(Y.seg[, i] ~ -1 + as.matrix(X.seg))
# allresids <- fm1$residuals
# allresids[abs(allresids)>1] <- NA
resids[, i] <- fm1$residuals
Bhat[, i] <- fm1$coefficients
}
Bhat
resids
Sigmahat <- cov(resids)
par_vals <- list(Sigmahat = Sigmahat, Bhat = Bhat)
result <- mapply(function(Y.seg,X.seg){
if (nrow(Y.seg) > (365)) {
# compare means of Y over the years
doit <- function() {
Bhat <- matrix(0, k, d)
resids <- matrix(0, n, d)
for (i in 1:d) {
fm1 <- lm(Y.seg[, i] ~ -1 + as.matrix(X.seg))
# allresids <- fm1$residuals
# allresids[abs(allresids)>1] <- NA
resids[, i] <- fm1$residuals
Bhat[, i] <- fm1$coefficients
}
Sigmahat <- cov(resids)
par_vals <- list(Sigmahat = Sigmahat, Bhat = Bhat)
return(par_vals)
}
par_vals_try <- try(
doit(),
silent = T
)
if (inherits(par_vals_try, "try-error")) {
par_vals <- list(Sigmahat = matrix(NA, d, d), Bhat = matrix(NA, k, d))
} else {
if (any(is.na(c(unlist(par_vals_try))))) {
par_vals <- list(Sigmahat = matrix(NA, d, d), Bhat = matrix(NA, k, d))
}
par_vals <- par_vals_try
}
} else {
par_vals <- list(Sigmahat = matrix(NA, d, d), Bhat = matrix(NA, k, d))
}
},Y.seg=Y.segs,X.seg=X.segs,SIMPLIFY = F)
result
result[[1]]
result[[2]]
result[[3]]
is.null(mvec)
60*24
