How can I find the 2 roots of a functions in R for different size of an argument?

[homerjaysimpson]

I am currently reading a statistical procedure that tests a null hypothesis. The site that I am reading is this one .

I have constructed a function in R that calculates the likelihood accordingly:

f = function(x,T=obs)  {
  X_a =  qchisq(0.05, df=1, ncp = 0, lower.tail = FALSE, log.p = FALSE) 
  p = 0.008
  2*log( (((T-x)/((1-p)*T))^(T-x))*( (x/(p*T))^x )  ) - X_a
}

Now the author on the site states that " Consider the example we looked at with the recommended standard coverage test. We implement a one-day 95% value-at-risk measure and plan to backtest it at the .05 significance level after 500 trading days, so q = 0.95 and α + 1 = 500. We calculate the ε = .05 quantile of the χ^2(1,0) distribution as 3.841. Setting this equal to [14.7](i.e the likelihood), we solve for X. There are two solutions: 16.05 and 35.11. Rounding down and up, respectively, we set x1 = 16 and x2 = 36. We will reject the value-at-risk measure if X ∉ [16, 36]."

one approach is that I can take :

K = function(obs){
f = function(x,T=obs)  {
  X_a =  qchisq(0.05, df=1, ncp = 0, lower.tail = FALSE, log.p = FALSE) 
  p = 0.008
  2*log( (((T-x)/((1-p)*T))^(T-x))*( (x/(p*T))^x )  ) - X_a
}
if (obs>=1000) {
  xlim.up = as.numeric(substr(obs,1,2))
  xlim.down = as.numeric(substr(obs,1,2))
} else {
  xlim.up = as.numeric(substr(obs,1,1))
  xlim.down = as.numeric(substr(obs,1,1))
}

x1 = uniroot(f, lower = 0,upper = xlim.up)  ### our 1st root 
x2 = uniroot(f, lower = xlim.down,upper =  80)  ### our 2sd root 
a =ceiling(x1$root)
b = floor(x2$root)
return(c(a,b))
}
K(2000)

with result:

K(2000)
>[1]  9 24

but when I lower the observation value to say 50 R reports me :

K(50)
 Error in uniroot(f, lower = xlim.down, upper = 80) : 
f() values at end points not of opposite sign

and it doesn't work.

My questions are:

1. how I can solve this equation in R and do the rounding?
2. For different values of observations starting from 50 up to 5000 how the function can be modified appropriately in order to report me the lower and upper bound as is being show in the table 14.4 in the site ?

Any help ?

[AlexisW]

First, note that if your question is whether a given observation is within bounds, you don't need to compute the interval. You can just do as the {segMGarch} package and compare your value of -2log(LAMBDA) to the qchisq.

Your code is mostly correct, except for the part where you try to take the first digits of obs as a breakpoint, there is little reason for it to work.

The problem with uniroot() is that it typically only works for a single root.

You can try with rootSolve::uniroot.all(). Make sure you check that you got 2 roots. For example:

f = function(x, T=obs, q = .95)  {
  X_a =  qchisq(0.05, df=1, ncp = 0, lower.tail = FALSE, log.p = FALSE) 
  p = 1 - q
  2*log( (((T-x)/((1-p)*T))^(T-x))*( (x/(p*T))^x )  ) - X_a
}

find_roots <- function(obs, q, upper_limit = 50){
  res <- rootSolve::uniroot.all(f, interval = c(0,upper_limit), n = 2, T = obs, q = q)
  if(length(res) != 2){
    stop("Failed to find 2 roots")
  }
  res
}

find_roots(obs = 500, q = 0.95)
#> [1] 16.05049 35.10622

^{Created on 2022-10-04 by the reprex package (v2.0.1)}

This reproduces the result from your book. However, the upper_limit parameter is critical, and this function will fail if you don't change this parameter depending on obs.

One "easy" way is to manually run it for a set of obs between 50 and 5000, and try to find a rule of thumb.

A more complex approach is to use the fact that you know the function shape. So you could write a function that does:

• start with a vector x
• compute f(x) and diff(f(x))
• ensure that the sign of diff(f(x)) is negative at the beginning, and positive at the end
  • if it's not negative at the beginning, take a more dense x
  • if it's not positive at the end, take a wider interval
• find where the sign of diff(f(x)) inverts
• call uniroot() on each of the segments

So a not-well-tested implementation can look like:

f = function(x, T=obs, q = .95)  {
  X_a =  qchisq(0.05, df=1, ncp = 0, lower.tail = FALSE, log.p = FALSE) 
  p = 1 - q
  2*log( (((T-x)/((1-p)*T))^(T-x))*( (x/(p*T))^x )  ) - X_a
}


find_two_roots <- function(FUN, density = .1, upper_limit = 50, ...){
  x <- density*(0:(upper_limit/density))
  diff_f <- diff(FUN(x, ...))
  
  if(diff_f[2] >= 0){
    #we miss the beginning
    find_two_roots(FUN, density = density*.1, upper_limit = upper_limit, ...)
  } else if(diff_f[length(diff_f)] <= 0){
    #we miss the end
    find_two_roots(FUN, density = density, upper_limit = upper_limit*10, ...)
  } else{
    # we have the correct beginning and end
    mid <- which(diff(sign(diff_f)) > 0)
    
    if(length(mid) != 1){
      stop("Incorrect inversions of the derivative: ", mid)
    }
    
    root1 <- uniroot(f, lower = 0, upper = x[mid], ...)$root
    root2 <- uniroot(f, lower = x[mid], upper = upper_limit, ...)$root
    return(c(floor(root1), ceiling(root2)))
  }
}
find_two_roots(f, T = 500, q = .95)
#> [1] 16 36

find_two_roots(f, T = 5000, q = .93)
#> [1] 315 386

^{Created on 2022-10-05 by the reprex package (v2.0.1)}

[homerjaysimpson]

Thank you for your answer. I see the problem is the interval.
I am thinking of the optimise function to directly search for a local minima or maxima instead of take the derivative of f and the then find the root.
In my real world problem q = .992.

ceiling(optimise(f,c(0,10),T = 50,q =.992,maximum =FALSE)$minimum)
floor(optimise(f,c(0,40),T = 50,q =.992,maximum =TRUE)$maximum)

if i use your solution with upper_limit 80 and for 500 and 5000 observations works but for less than 500 it doesn't

> find_two_roots(f, T = 80, q = .992)
 Error in uniroot(f, lower = 0, upper = x[mid], ...) : 
f() values at end points not of opposite sign 
> find_two_roots(f, T = 500, q = .992)
[1] 0 9
> find_two_roots(f, T = 5000, q = .992)
[1] 28 53

and help ?