Function to output forecasting skill. Introduces new type of forecasting skill. Includes unit tests.

fskill.r

mse = function(v1,v2){
  se = (v2 - v1)^2
  mse = mean(se)
  return(mse)
}

# See http://en.wikipedia.org/wiki/Forecast_skill for info on forecasting skill. 
# One thing I don't like about the vanilla method is that the domain is (-Inf,1] which is kind of hard to understand.
# The Rossignol method, or Rossignol's Modified Forecasting Skill (RMFS) changes this to be [-1,1] which is easier to understand. Use it by setting type="rossignol".
fskill = function(data,type="vanilla"){
  # check inputs
  if(ncol(data) < 3){
    print("Format of data must be data frame with at least 3 columns. First column are the actual values, second column is base predictions, other columns are new predictions.")
    stop()
  }
  
  if(any(apply(data,2,class) != "numeric")){
    print("One column in input contains non-numeric values.")
    stop()
  }
  
  if(!(type %in% c("vanilla","rossignol"))){
    print(paste("Invalid type:",type,"using default 'vanilla'."))
  }
  
  # main body
  mse.base = mse(data[,1],data[,2])
  ret.vals = c()
  for(i in 3:ncol(data)){
    mse.new = mse(data[,1],data[,i])
    if(type == "vanilla"){
      skill = 1 - mse.new/mse.base
    }
    if(type == "rossignol"){
      if(mse.new == mse.base) skill = 0
      if(mse.new < mse.base) skill = 1 - mse.new/mse.base
      if(mse.new > mse.base) skill = mse.base/mse.new - 1
    }
    ret.vals = c(ret.vals,skill)
  }
  return(ret.vals)
}

require(testthat)
set.seed(1)
data = data.frame(actual = runif(100),base = runif(100),new1 = runif(100),new2 = runif(100))
test_that("vanilla fskill returns correct values",expect_equal(round(fskill(data),2),c(0.06,-0.27)))
data = data.frame(actual = runif(100),base = runif(100),new1 = runif(100,1,2),new2 = runif(100))
test_that("rossignol fskill returns correct values",expect_equal(round(fskill(data,type="rossignol"),2),c(-0.86,-0.03)))