Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R Fit linear line or exponential curve through peaks and lows of timeseries

Tags:

plot

r

line

I would like to plot a linear or log line through the peaks and lows in a times series with decreased peaks and lows like in the images. How do I do that?

In oher words: a line with its origin at the data series' peak, with the most negative slope possible, that doesn't intersect the data series at any other point. Which also accounts for the lows. Preferably both a lineair line as an exponential line.

The red lines I had drawn manually cause I dont know how to get them via R, hence thats the goal of the question. The following code gives the timeseries.

set.seed(123)

# Generate time series 
n <- 10000
x <- cumsum(rnorm(n, 0, 1))
sigma <- seq(2, 0.1, length.out = n)
x <- x * sigma

# Plot the time series
plot(x, type = "l")

enter image description here enter image description here

like image 848
H. berg Avatar asked Nov 17 '25 11:11

H. berg


2 Answers

1) max/min slope Determine the maximum slope between the maximum and points to the right of the maximum and the minimum slope from the minimum to points to the right of the minimum.

plot(x, type = "l")
n <- length(x)

w.mn <- which.min(x)
x1 <- x[w.mn:n]
m1 <- min((x1[-1] - min(x1)) / ((w.mn+1):n - (w.mn)))
segments(w.mn, min(x), n, min(x) + m1 * (n - w.mn), col = "red", lwd = 2)

w.mx <- which.max(x)
x2 <- x[w.mx:n]
m2 <- max((x2[-1] - max(x2)) / ((w.mx+1):n - (w.mx)))
segments(w.mx, max(x), n, max(x) + m2 * (n - w.mx), col = "red", lwd = 2)

screenshot

2) quantile regression. Another way to do this is to use quantile regression. We show the lines in red and the log upper curve in green. The lower curve does not really make sense for log due to negative values in x.

library(quantreg)

xx <- seq_along(x)
plot(x ~ xx, type = "l")

fm <- rq(x ~ xx, tau = 0:1)
yy <- predict(fm)
lines(yy[, 1] ~ xx, col = "red", lwd = 2)
lines(yy[, 2] ~ xx, col = "red", lwd = 2)

fm.log <- rq(log(x) ~ xx, tau = 1)
yy.log <- cbind(1, xx) %*% coef(fm.log)
lines(exp(yy.log) ~ xx, col = "green", lwd = 2)

screenshot

like image 189
G. Grothendieck Avatar answered Nov 20 '25 00:11

G. Grothendieck


Based on @zephyrl's suggestion, you could do a grid search over the slopes that have a starting point at the peak (or valley). The stopping criterion would be intersection of any point in the series. Here's a function that would do it. The series argument is the series of data points, the time argument is for times associated with those series (to make the plot, step is the increment by which the slope is changed in each iteration and log is logical indicating whether the time variable should be logged.

set.seed(123)

# Generate time series 
n <- 10000
x <- cumsum(rnorm(n, 0, 1))
sigma <- seq(2, 0.1, length.out = n)
x <- x * sigma
find_lines <- function(series, time=NULL, step=1e-5, tol=1e-4, log=FALSE, ...){
  if(is.null(time)){
    time <- seq_along(series)
  }
  tmp <- data.frame(time = time, x=x)
  mx_tmp <- tmp[which.max(tmp$x), ]
  mx_tmp$x <- mx_tmp$x + tol
  mn_tmp <- tmp[which.min(tmp$x), ]
  mn_tmp$x <- mn_tmp$x - tol
  b_top <- 0
  b_bot <- 0
  if(!log){
    a_top <- mx_tmp$x - b_top*mx_tmp$time 
    a_bot <- mn_tmp$x - b_bot*mn_tmp$time 
    intersect_top <- any(tmp$x > b_top*tmp$time + a_top)    
    intersect_bot <- any(tmp$x < b_bot*tmp$time + a_bot)    
    while(!intersect_top){
      b_top <- b_top - step    
      a_top <- mx_tmp$x - b_top*mx_tmp$time 
      intersect_top <- any(tmp$x > b_top*tmp$time + a_top)    
    }  
    b_top <- b_top + step
    a_top <- mx_tmp$x - b_top*mx_tmp$time 
    while(!intersect_bot){
      b_bot <- b_bot + step    
      a_bot <- mn_tmp$x - b_bot*mn_tmp$time 
      intersect_bot <- any(tmp$x < b_bot*tmp$time + a_bot)    
    }  
    b_bot <- b_bot - step
    a_bot <- mn_tmp$x - b_bot*mn_tmp$time 
  }else{
    a_top <- mx_tmp$x - b_top*log(mx_tmp$time)
    a_bot <- mn_tmp$x - b_bot*mn_tmp$time 
    intersect_top <- any(tmp$x > b_top*log(tmp$time) + a_top)    
    intersect_bot <- any(tmp$x < b_bot*log(tmp$time) + a_bot)    
    while(!intersect_top){
      b_top <- b_top - step    
      a_top <- mx_tmp$x - b_top*log(mx_tmp$time)
      intersect_top <- any(tmp$x > b_top*log(tmp$time) + a_top)    
    }  
    b_top <- b_top + step
    a_top <- mx_tmp$x - b_top*log(mx_tmp$time)
    while(!intersect_bot){
      b_bot <- b_bot + step    
      a_bot <- mn_tmp$x - b_bot*log(mn_tmp$time) 
      intersect_bot <- any(tmp$x < b_bot*log(tmp$time) + a_bot)    
    }  
    b_bot <- b_bot - step
    a_bot <- mn_tmp$x - b_bot*log(mn_tmp$time)
  }
  return(list(bottom = list(a = a_bot, b=b_bot), 
              top = list(a = a_top, b=b_top), 
              log=log))
  
}
f <- find_lines(x)
plot(x, type="l")
abline(f$top$a, f$top$b, col="red")
abline(f$bot$a, f$bot$b, col="red")
f <- find_lines(x, log=TRUE, step=1e-1)
lines(1:length(x), f$top$a + f$top$b*log(1:length(x)), col="green")
lines(1:length(x), f$bot$a + f$bot$b*log(1:length(x)), col="green")

Created on 2023-03-30 with reprex v2.0.2

like image 44
DaveArmstrong Avatar answered Nov 20 '25 01:11

DaveArmstrong