.Rhistory

return(prod(product_runs))
}
old_new_product_death = function(k){
for (j in k){
product_runs = append(product_runs,(n-3*ma-(q1-1)*(2*ma)+j))
}
return(prod(product_runs))
}
prior = NA
birth_old_to_new_ratio = (factorial(num_of_bkpts+1)*q2(dpois(k_ends_new,lambda)))/(old_new_product_birth(num_of_bkpts)*prior)
birth_new_to_old_ratio = (factorial(num_of_bkpts)*q1(dpois(k_ends,lambda)))/(product_death(num_of_bkpts)*prior)
death_new_to_old_ratio = (factorial(num_of_bkpts)*q2*(dpois(k_ends,lambda)))/(product_death(num_of_bkpts)*prior)
death_old_to_new_ratio = (factorial(num_of_bkpts-1)*make_k*(dpois(k_ends_new,lambda)))/(old_new_product_death(num_of_bkpts)*prior)
#
#end of birth and death ratios
delta_bic = (-2*new_loglik + log(n)*(length(k_ends_new)-1)*(3+ma)) - (-2*old_loglik + log(n)*(length(k_ends)-1)*(3+ma))
ratio = (-1*delta_bic/2) + (log(q1*dnorm(length(k_ends_new)-2,epsilon_list[ma],init_sigma2)) - log(q2*dnorm(length(k_ends)-2,epsilon_list[ma],init_sigma2))) #both dnorms are the same, make sure that this is correct
u_ratio = log(runif(1)) #random number from 0 to 1 taken from a uniform distribution and then log transformed
if(abs(delta_bic) == Inf){ #safe guard against random models creating infinite ratios
k_ends <<- k_ends #old
} else if(ratio > u_ratio) {
k_ends <<- k_ends_new #new
} else {
k_ends <<- k_ends #old
}
if(progress == TRUE){
setTxtProgressBar(burn_progress, i)
}
}
#initializing matrices/storage objects for final Metropolis-Hasting
all_k_best = data.frame(matrix(ncol=(length(k_ends)-2),nrow=0))
if(fit_storage == TRUE){
bar_v = 0
bar_beta = 0
fit = 0
all_fits = data.frame()
all_MSE = data.frame()
}
all_BIC = data.frame()
accept_count = 0
#setting up counters for final Metropolis-Hasting
type = "0"
a.count <<- 0
s.count <<- 0
m.count <<- 0
j.count <<- 0
add.accept.count <<- 0
sub.accept.count <<- 0
move.accept.count <<- 0
jiggle.accept.count <<- 0
#getting constants for qs for final Metropolis-Hasting
starting_bkpts = length(k_ends) - 1 #most probable number of breakpoints based on starting info
starting_nfree = length(freeObservations(k_ends,ma))
burn_in = 50
make_murder_p = 0.5
percent = 0.02
lambda = 1
tao = abs(1)
bs(1), alpha = 0, mu = 1.5, jump_p = 0.25, ma = 1, progress = TRUE, fit_storage = TRUE)
jump_p = 0.25
ma = 1
mu = 1.5,
mu = 1.5
alpha = 0
fit_storage = TRUE
progress = TRUE
burn_in = 50
make_murder_p = 0.5
lambda = 1
percent = 0.02
y = rnorm(1,0,1)
mu = 1
alpha = .5
epsilon = c(rnorm(1,0,1))
epsilon = function(t){ # t MUST be greater than or equal to 2
for (i in 2:t){
epsilon[i-1] = rnorm(1,0,1)
}
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon(2) # + .............
}
y[i]= mu + epsilon[1] + alpha*epsilon(2) # + .............
epsilon[1]
epsilon = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 2:t){
epsilon_t[i-1] = rnorm(1,0,1)
}
}
y[i]= mu + epsilon + alpha*epsilon(2) # + .............
epsilon = rnorm(1,0,1)
y[i]= mu + epsilon + alpha*epsilon(2)
epsilon = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 2:t){
epsilon_t[i-1] = rnorm(1,0,1)
}
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon(2) # + .............
}
epsilon = rnorm(1,0,1)
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon(2) # + .............
}
epsilon = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 2:t){
epsilon_t[i-1] = rnorm(1,0,1)
}
}
mu = 1
y = rnorm(1,0,1)
alpha = .5
epsilon = rnorm(1,0,1)
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon(2) # + .............
}
mu = 1
y = rnorm(1,0,1)
alpha = .5
epsilon = rnorm(1,0,1)
epsilon_t = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 2:t){
epsilon_t[i-1] = rnorm(1,0,1)
}
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
mu = 1
y = rnorm(1,0,1)
alpha = .5
epsilon = rnorm(1,0,1)
epsilon_t = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 1:t){
epsilon_t[i] = rnorm(1,0,1)
}
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
mu = 1
y = rnorm(1,0,1)
alpha = .5
epsilon = rnorm(1,0,1)
for (i in 1:t){
epsilon_t[i] = rnorm(1,0,1)
}
epsilon_t = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 1:t){
epsilon_t[i] = rnorm(1,0,1)
}
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
epsilon_t(2)
epsilon_t(2)
mu
epsilon
alpha
epsilon_t = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 1:t){
epsilon_t[i] = rnorm(1,0,1)
}
return(epsilon_t[t])
}
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
epsilon_t[i] = rnorm(1,0,1)
epsilon_t = function(t){ # t MUST be greater than or equal to 2
epsilon_t = c(epsilon)
for (i in 1:t){
epsilon_t[i] = rnorm(1,0,1)
}
return(epsilon_t[t])
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
beta=-1
y=rnorm(1,0,1)
for(i in 2:32){y[i]=beta*y[i-1]+rnorm(1,0,1)}
for(i in 33:90){y[i]=beta*y[i-1]+rnorm(1,0,3)}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
for (i in 1:90){
y[i]= mu + epsilon + alpha*epsilon_t(2) # + .............
}
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  ts() #time series
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
splice_data = function(current_data){ #Why does this function need to be called for each bkpt?
for (j in 1:ncol(current_data)){
if (j == 1){ #if j is time
append(new_time,c(current_data[,j][1:i]))
} else if (j==2) {
append(new_data,c(current_data[,j][1:i]))
}
}
new_start = data.frame(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
return(new_start)
}
spliced_data = splice_data(current_data)
#append to list of arima'ed time series
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
}
return(new_fits)
}
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
library(devtools)
#install_version("FitAR", version = "1.94", repos = "http://cran.us.r-project.org")
#library(FitAR)
library(forecast)
library(strucchange)
library(stats)
library(MASS)
par(mar=c(1,1,1,1))
setwd("/Users/mellm/github/REU2022/test_Cases")
suicide<-read.csv("suicideRatesData.csv")
attach(suicide)
setwd("/Users/mellm/github/REU2022/")
source("ar_bai_perron.r")
source("data_for_trials.r")
source("Variations_of_bars/fixedbaar1.r")
#fit data with ARIMA
single_model <- arima(suicide$X15.24.years, order=c(2,0,0))
single_fitted <- fitted(single_model)
test_data_45 = function(){
data_45=suicide$X15.24.years
time = c(1:37)
test_data_45 = data.frame(time, data_45)
return(test_data_45)
}
iterations=1000
runs=10
L=matrix(NA,nrow=iterations,ncol=runs)
A=NA
M=NA
B=NA
for(i in 1:runs){
current_data = test_data_45()
break_p = breakpoints(current_data[,2] ~ current_data[,1], breaks = 3, h = 0.1)
starting_breakpoints = break_p$breakpoints
test1=baar(starting_breakpoints,1:37,current_data[,2],iterations, 15, jump=0.25, ar=3)
print(i)
L[,i]=test1$NumBkpts
M[i]=length(starting_breakpoints[!is.na(starting_breakpoints)])
A=c(A,(unlist(test1$Breakpoints))) #Baar location of breakpoints
B=c(B,(starting_breakpoints)) #B-P
}
A2=A[!is.na(A)]
B2=B[!is.na(B)]
mean(L)
sd(L)
mean(M)
sd(M)
current_data = ts(current_data)
A3 = sort(unique(A2))
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
new_start = data.frame(new_time,new_data)
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  ts() #time series
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
splice_data = function(current_data){ #Why does this function need to be called for each bkpt?
for (j in 1:ncol(current_data)){
if (j == 1){ #if j is time
append(new_time,c(current_data[,j][1:i]))
} else if (j==2) {
append(new_data,c(current_data[,j][1:i]))
}
}
new_start = data.frame(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
return(new_start)
}
spliced_data = splice_data(current_data)
#append to list of arima'ed time series
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
}
return(new_fits)
}
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
current_data
current_data = current_data
current_data = test_data_45()
current_data = current_data
as.data.frame(current_data)
is.data.frame(current_data)
is.data.frame(single_model)
is.ts(single_model)
typeof(single_model)
current_data = as.list(current_data)
current_data
current_data[,1]
current_data[1]
current_data[1][1:5]
current_data[1][2]
current_data[1][[2]]
current_data[[1]]
current_data[[1]][1:5]
is.list(current_data[[1]])
typeof(current_data[[1]])
c(current_data[[1]][1:5])
typeof(current_data[[1]][1:5])
typeof(c(current_data[[1]][1:5]))
as.list(current_data[[1]][1:5])
append(current_data[[1]][1:5], 3)
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  list() #list
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
for (j in 1:ncol(current_data)){
if (j == 1){ #if j is time
append(new_time,current_data[[j]][1:i])
} else if (j==2) {
append(new_data,current_data[[j]][1:i])
}
}
}
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
}
return(new_fits)
new_fits = c() #list of new fitted values
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  list() #list
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
for (j in 1:ncol(current_data)){
if (j == 1){ #if j is time
append(new_time,current_data[[j]][1:i])
} else if (j==2) {
append(new_data,current_data[[j]][1:i])
}
}
}
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
}
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  list() #list
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
for (j in 1:ncol(current_data)){
if (j == 1){ #if j is time
append(new_time,current_data[[j]][1:i])
} else if (j==2) {
append(new_data,current_data[[j]][1:i])
}
}
}
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
return(new_fits)
}
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
ncol(current_data)
length(current_data)
new_baar_fit = function(current_data,breakpoints){
new_fits = c() #list of new fitted values
spliced_fitted_arima =  list() #list
new_time = c() #list of time vectors for each bkpt (for spliced_data)
new_data = c() #list of data_45 vectors for each bkpt (for spliced_data)
for (i in breakpoints){
for (j in 1:length(current_data)){
if (j == 1){ #if j is time
append(new_time,current_data[[j]][1:i])
} else if (j==2) {
append(new_data,current_data[[j]][1:i])
}
}
}
new_start = list(new_time,new_data) #data frame of time and data_45 vectors for each bkpt
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here
new_fits = fitted(spliced_fitted_arima)
return(new_fits)
}
#WORK IN PROGRESS ---------------------------------------------------
fits = new_baar_fit(current_data,A3)
new_start = cbind(c(1:length(as.numeric(new_time))), as.numeric(new_data)) #combining time and data inputs
View(new_start)
new_time = c(1,2,4) #list of time vectors for each bkpt (for spliced_data)
new_data = c(1,2,3) #list of data_45 vectors for each bkpt (for spliced_data)
new_start = cbind(c(1:length(as.numeric(new_time))), as.numeric(new_data)) #combining time and data inputs
new_start = cbind(c(as.numeric(new_time)), as.numeric(new_data)) #combining time and data inputs
new_start = ts(c(as.numeric(new_time)), as.numeric(new_data)) #combining time and data inputs
spliced_fitted_arima = arima(new_start$new_data,order = c(2,0,0)) #Something wrong here