[R-es] Prueba de Thom y series temporales homogéneas de prueba

[Rubén Gómez Antolí]

Hola a todos:

Hace unos meses estuve preguntando sobre la prueba de Thom para 
comprobar la homogeneidad de series temporales sobre datos de radiación 
e insolación (después he visto que también se usa con otras variables 
meteorológicas) [0].

El caso es que por distintas razones no había podido avanzar en el 
asunto y ahora vuelvo sobre él.

Creé una función para realizar la prueba de Thom pero el caso es que 
series temporales de radiación que «deberían» pasar la prueba no lo 
hacen. (Digo deberían porque he encontrado múltiples artículos donde 
mencionan los mismos datos que yo uso (AEMET) y pasan la prueba).

Como de alguna forma quiero comprobar que la función esta bien 
construida he intentado buscar alguna documentación para crear series 
temporales homogéneas pero no he conseguido dar con nada para hacerlo 
(puede deberse también a mi torpeza, claro).

Expongo el código a continuación:

--------------- código R -------------
contar.pasos<-function(vector,valor.cruce) {
     # Función que cuenta los cruces de un vector por un
     # valor de referencia (valor.cruce)

     paso<-0
     i<-1

     while ( i <= (length(vector)-1) ) {

       cruce.corto<-(vector[i] < valor.cruce & vector[(i+1)] > 
valor.cruce ) |
	  (vector[i] > valor.cruce & vector[(i+1)] < valor.cruce )

       if ( vector[i] == valor.cruce & i!=1 ) {
       # La función esta en el valor de cruce y no es el primer valor,
       # se presentan varios casos

	baja<-( (vector[i-1] > vector[i]) & (vector[i] > vector[i+1]) )
	sube<-( (vector[i-1] < vector[i]) & (vector[i] < vector[i+1]) )

	if (sube | baja ) {pcaso<-TRUE} else {pcaso<-FALSE}
	# La función cruza el valor de referencia.

	if (vector[(i-1)] == vector[i]) {
	  # El valor anterior de la función también esta en el valor de cruce,
	  # entonces comprueba como era la función después del tramo horizontal
	  # y como evoluciona después
	  z<-i-1
	  while (vector[z] == vector[i] & z>1) {z<-z-1}
	  # Retrocedemos hasta el último valor distinto sin pasar
	  # del primer valor de la función.
	  menor<-FALSE
	  mayor<-FALSE
	  if (vector[z] < vector[i]) { menor<-TRUE } else {mayor<-TRUE}
	  if (
	    ( menor & (vector[i] < vector[(i+1)]) ) |
	    ( mayor & (vector[i] > vector[(i+1)]) )
	  ) {scaso<-TRUE}
	} else {scaso<-FALSE}
	
	cruce.largo<-(pcaso | scaso)

       } else {cruce.largo<-FALSE}

       excepcion<-if (
	i==1 | (vector[i] == valor.cruce & !cruce.largo)
	) {TRUE} else {FALSE}
	
       if ( (cruce.corto | cruce.largo) & !excepcion ) {paso<-paso+1}
	# Si la función cruza y no hay excepción

       i<-i+1
     }

     return(paso)
}

pdThom<-function(vector,na.rm=F) {
     # Prueba de Thom, de las rachas o de las alternancias. (1.966)
     # Prueba de homogeneidad absoluta aconsejada para distribuciones no
     # conocidas o no normales.
     # Depende de la función contar.pasos

     longitud<-length(vector)
     mediana<-median(vector,na.rm=na.rm)
     rachas<-contar.pasos(vector,mediana)

     resultado<-(rachas-((longitud+2)/2))/
       sqrt((longitud*(longitud-2))/(4*(longitud-1)))

     return(abs(resultado))
}

library(solaR)

SIAR <- read.csv('http://solar.r-forge.r-project.org/data/SIAR.csv')
#
Estaciones.Almeria<-subset(SIAR,Provincia == "Almería")
#
# Datos estación de La Mojonera
# Almería, provincia = 4, estación = 1
Estacion.La.Mojonera<-readSIAR(4,1,"01/01/1970","31/12/2011")
#
# Datos estación de Adra
# Estación = 10
Estacion.Adra<-readSIAR(4,10,"01/01/1970","31/12/2011")
#
# Datos estación Almería
# Estación = 2
Estacion.Almeria<-readSIAR(4,2,"01/01/1970","31/12/2011")
#
# Radiación diaria en Adra.
Adra.radiacion.diaria<-getG0(Estacion.Adra)

# Radiación diaria en La Mojonera
Mojonera.radiacion.diaria<-getG0(Estacion.La.Mojonera)

# Radiación diaria en Almería
Almeria.radiacion.diaria<-getG0(Estacion.Almeria)

# Temperatura diaria en Almeria
Almeria.temp<-Estacion.Almeria en data$TempMedia

## Pruebas de homogeneidad
which(is.na(as.vector(Almeria.temp)))
# Resultado: integer(0), no hay valores perdidos

pdThom(as.vector(Almeria.temp))
# Resultado: [1] 59.20887, bastante alejado de 2,58

which(is.na(as.vector(Adra.radiacion.diaria)))
#  [1]    1   11  825 1435 1449 1450 1451 1452 1453 1454 1756 1868 1895 
# 1896 1897
# [16] 1903 1947 1948 1949 1950 1951 1952 1954 3465 3614

pdThom(as.vector(Adra.radiacion.diaria)[12:824])
# [1] 22.63508

which(is.na(as.vector(Almeria.radiacion.diaria)))
# [1]  886 1838 3221 3222 3223 3224 3227 3402 3470 3592 3693

pdThom(as.vector(Almeria.radiacion.diaria)[1839:3220])
# [1] 28.201

library(TSA)
data(tempdub)
pdThom(tempdub)
# [1] 8.195372
data(larain)
pdThom(larain)
# [1] 0.4683109
------------------- Fin del código R -----------

Como veis la única serie que «pasa» la prueba de Thom es la de lluvia en 
Los Angeles del paquete TSA, pero estoy desorientado, no me lo creo.

¿Que opináis?

Gracias por adelantado.

Un saludo.

R. Gómez

[0] https://stat.ethz.ch/pipermail/r-help-es/2012-February/003236.html
-- 
Libertad es poder elegir en cualquier momento. Ahora yo elijo GNU/Linux,
para no atar mis manos con las cadenas del soft propietario.
---------
Desde El Ejido, en Almería, usuario registrado Linux #294013
http://www.counter.li.org