#Transition probability/ Overgangssannsynlighet
#Fix G and R, my.kol(x,t) is the force of mortality for an x-year old
#Norwegian at calendar year t (bigger than 2013)

source("K2013.r")

x0 <- 35 #age at year Y
Y <- 2023 #current year of choice
G <- 0 #gender
R <- 0 #risk

f <- function(s) my.kol(x0+s,Y+s)
h <- 0.01 #step size of partition

Riemann.left <- function(a,b){
  ti <- seq(a,b,by=h)
  N <- length(ti)-1
  fi <- sapply(ti,f)
  return(h*sum(fi[1:(N-1)]))
}

Riemann.mid <- function(a,b){
  ti <- seq(a,b,by=h)
  N <- length(ti)-1
  fi <- sapply(ti+(h/2),f)
  return(h*sum(fi[1:(N-1)]))
}

Riemann.right <- function(a,b){
  ti <- seq(a,b,by=h)
  N <- length(ti)-1
  fi <- sapply(ti,f)
  return(h*sum(fi[2:N]))
}

trapezoidal <- function(a,b){
  ti <- seq(a,b,by=h)
  N <- length(ti)-1
  fi <- sapply(ti,f)
  return((h/2)*sum(fi[1:(N-1)] + fi[2:N]))
}

Simpson <- function(a,b){
  ti <- seq(a,b,by=h)
  N <- length(ti)-1
  fi <- sapply(ti,f)
  aux1 <- aux2 <- 0
  for(i in 2:(N-1)){
    if( i%%2 == 0 ){ aux1 <- aux1+2*f(ti[i]) }
    if( i%%2 == 1 ){ aux2 <- aux2+4*f(ti[i]) }
  }
  return((h/3)*(fi[1]+aux1+aux2+fi[N]))
}

#Compute survival probability with you favourite method
p.K13.riemann.left <- function(n,m) exp(-Riemann.left(n,m))
p.K13.riemann.mid <- function(n,m) exp(-Riemann.mid(n,m))
p.K13.riemann.right <- function(n,m) exp(-Riemann.right(n,m))
p.K13.trap <- function(n,m) exp(-trapezoidal(n,m))
p.K13.simp <- function(n,m) exp(-Simpson(n,m))
p.K13.Rintegrate <- function(n,m) exp(-integrate(f,n,m)$value)
library(cubature)
p.K13.Rcubature <- function(n,m) exp(-hcubature(f, n, m)$integral)

G <- 0
R <- 0
x0 <- 35
n <- 0
m <- 60
p.K13.riemann.left(n,m)
p.K13.riemann.mid(n,m)
p.K13.riemann.right(n,m)
p.K13.trap(n,m)
p.K13.simp(n,m)
p.K13.Rintegrate(n,m)
p.K13.Rcubature(n,m)

G <- 1
R <- 0
x0 <- 35
n <- 0
m <- 60
p.K13.riemann.left(n,m)
p.K13.riemann.mid(n,m)
p.K13.riemann.right(n,m)
p.K13.trap(n,m)
p.K13.simp(n,m)
p.K13.Rintegrate(n,m)
p.K13.Rcubature(n,m)