fce <- Vectorize(function(t){
  integrand <- function(x) 2*(x-t)/(1+(x-t)^2)*dcauchy(x);
  integrate(integrand, -50, 50)$value 
})

t <- seq(-10,10, by=0.1);
plot(t, fce(t), type="l")

x <- rcauchy(5);
psi.stat <- Vectorize(function(t) sum((x-t)/(1+(x-t)^2)));

t <- seq(-15,15, by=0.1);
plot(t, psi.stat(t), type="l")
abline(h=0, lty="dotted")

# diff(psi.stat)
# 
# [1] -0.5667116 -0.3193519 12.5888410  0.6558730  0.9608010

# # 
# [1]   1.5210493 -11.8565599  -3.5805595   0.1394367  -0.9303372


DF <- 3:10; 

plot(DF, 1/(4*dt(0, df=DF)^2), type="l", ylim=c(0,2)); lines(DF, DF/(DF-2))


sigma.fce <- Vectorize(function(DF){
  integrand <- function(x) x^2/(1+x^2)^2*dt(x,df=DF);
  integrate(integrand, -10, 10)$value 
})

gamma.fce <- Vectorize(function(DF){
  integrand <- function(x) (1/(1+x^2)-2*x^2/(1+x^2)^2)*dt(x,df=DF);
  integrate(integrand, -10, 10)$value 
})


lines(DF, sigma.fce(DF)/(gamma.fce(DF)^2), col="blue")


