Procesos Estoc´ asticos asticos Taller No 2: An´ alisis procesos ARMA con Matlab alisis Octubre 2007
1.
Funcio unciones nes Matl Matlab ab para para an´ analisis a ´lisis de procesos ARMA
1. Funciones unciones para identifica identificaci´ ci´on on de los ´ordenes ordenes p y q: a ) autoco autocorr: rr: calcu calcula la y grafic graficaa la funci funci´´on on de autocorrelaci´on on estimada. parrcorr: rr: calc calcula ula y grafi grafica ca la func funci´ i´on on de autocorrela autocorrelaci´ ci´on on parcial estimada. b ) parrco 2. Funciones unciones para estimaci´ estimaci´ on on y ajuste: a ) arma armaba bat: t: func funci´ i´on on para identificar la pareja (p,q) que produce el modelo con menor criterio de informaci´on on de Akaike (AIC). Ver pag. web de H. Hurd. 1 b ) arma armax: x: func funci´ i´on on que estima los par´emetros emetros del modelo ARMA(p,q). Matlab provee varias varias funciones para estimaci´on on de modelos pero solamente se considerar´a ´esta. ta . c ) resi resid: d: func funci´ i´ on para calcular los residuos del modelo con el fin de poder on realizar pruebas de hip´otesis otesis para determinar si es ruido blanco, como la prueba de Ljung-Box. Los residuos son estimaciones de Z n . d ) lbqte lbqtest st:: func funci´ i´ on para realizar la prueba de Ljung-Box. on e ) compar compare: e: funci´ funci´ on para examinar la calidad de los pron´osticos que se pueden on hacer con el modelo ajustado con el fin de determinar si ´este este es adecuado o n´o. o. Notaci´ on on Matlab. En lugar de la letra L, Matlab utiliza q−1 , luego q −1(X n ) = X n−1 . Por ejemplo, un modelo ARMA(4,2) se expresa en Matlab as´ as´ı: ı: A(q −1)X n = B (q −1)Z n 1
http://www.stat.unc.edu/faculty/hurd/stat185Data/progdoc.html
1
A(q−1) = 1 − 0,8875q−1 + 0,6884q −2 − 0,8956q −3 + 0,7022q−4 B(q−1) = 1 + 0,2434q−1 + 0,4649q −2 1.0.1.
Pasos para el an´ alisis con las funciones de Matlab
1. Graficar la trayectoria del proceso, la FAC estimada, la FACP estimada y el Variograma. Las instrucciones a continuaci´on muestran c´omo hacerlo. figure(1) subplot(2,2,1), plot(x); ylabel(’Xn’) title(’Trayectoria’) [fac_y,m]=autocorr(x,[],2); subplot(2,2,2), autocorr(x,[],2) title(’fac’); [facp_y, mp] = parcorr(x,[],2); subplot(2,2,3), parcorr(x,[],2) title(’facp’); v = (fac_y(1)-fac_y)/(fac_y(1)-fac_y(2)); subplot(2,2,4), stem(m,v); grid title(’Variograma’)
2. Despu´ es de una posible identificaci´on del tipo de proceso se procede a especificar los ´ordenes p y q del proceso. En caso de no ser posible identificar un AR ´o un MA, se toma inicialmente el rango p, q = 1, 2, 3, 4, 5, 6, y se corre la funci´on “armabat”. Antes de aplicar la funci´on es conveniente restar la media para obtener un proceso de media cero, asumiendo que el proceso es estacionario en covarianza: % eliminar la media xt = x -mean(x); % explora el orden pvec = [1 2 3 4 5 6]; qvec = [1 2 3 4 5 6]; [mbest,minaic,pbest,qbest]=armabat(xt,pvec,qvec); pbest qbest
2
En las variables “pbest” y “qbest” est´an los valores de los ´ordenes p y q que mejor describen el proceso. 3. Estimaci´on los par´ametros del modelo ARMA(p,q). Para esto se utiliza la funci´on “armax” con la pareja ( p,q) escogida en el punto anterior. Esta instrucci´on crea un objeto de nombre (por ejemplo)“armapq” que contiene varios campos con informaci´on sobre el modelo estimado. armapq = armax(xt,[pbest qbest]); present(armapq) armapq.a armapq.c
4. Pruebas de significaci´on de los Par´ametros. La instrucci´on tcrit = armapq.ParameterVector./sqrt(diag(armapq.CovarianceMatrix))
calcula un vector de cocientes t j = φ j /S φj , donde S φj es la desviaci´on est´andar
del coeficiente φ j . Cada valor t j es el valor de un estad´ıstico t de Student que sirve para probar la hip´otesis de que los coeficientes son significativamente diferentes de cero, es decir, se acepta la hip´otesis H o : φ j = 0 al nivel de 5% si observamos que |t j | < 1,96; en caso contrario se rechaza.
Tambi´en se puede calcular el valor estimado de la varianza del ruido blanco, σ 2, mediante la instrucci´on: armapq.NoiseVariance. 5. Para completar el an´alisis es necesario chequear si los residuos del modelo ajustado son ruido blanco. Los residuos son valores estimados del proceso Z n . La forma de hacerlo es calculando la fac y la fac parcial con los residuos. Si los residuos resultan ruido blanco ambas funciones deben mostrar todos los valores dentro de las bandas de Bartlett. El c´alculo de los residuos se puede hacer con los siguientes comandos. dato = iddata(xt); rarmapq = resid(armapq,dato); et = rarmapq.OutputData; [H,P,Qstat,CV] = lbqtest(et, [20 25]’, 0.05); [H,P,Qstat,CV] figure(3)
3
subplot(2,2,1), plot(et); title(’Residuos’) [fac_x,m] = autocorr(et,30,[],2); subplot(2,2,2), autocorr(et,30,[],2); title(’fac muestral’) subplot(2,2,3), parcorr(et,30,[],2); title(’facp muestral’) v = (fac_x(1)-fac_x)/(fac_x(1)-fac_x(2)); subplot(2,2,4), stem(m,v); grid title(’Variograma’)
6. Una manera de chequear si el modelo propuesto ajust´o bien los datos es ajustar el modelo con la primera mitad de los datos y utilizar la parte restante para comparar con los pron´osticos a un paso: se compara X n con el pron´ostico de X n realizado con el modelo. La funci´on “compare” de Matlab hace este c´alculo. figure(4) % uso de la funcion "compare" mitad = floor(length(xt)/2); ye = xt(1:mitad); yv = xt(mitad+1:end); model= armax(ye,[pbest qbest]); compare(yv,model,1);
Ejemplo 1.1. En la Figura (1) se puede ver la gr´ afica de la serie: entrada de gas
metano en un horno a gas ( Methane input into gas furnace: cu. ft/min. Sampling interval 9 seconds. Source: Box and Jenkins ( 2)), se aprecian la fac y la fac parcial estimadas. Seg´ un lo explicado acerca de la identificaci´ on de modelos tipo ARMA, corresponder´ıa a un proceso autorregresivo AR(3). La estimaci´ on con armabat y armax se realiz´ o de acuerdo a las indicaciones anteriores y se obtuvo un modelo ARMA(1,4)(!!!). El modelo ajustado para la serie X n − µ es: Discrete-time IDPOLY model: A(q)y(t) = C(q)e(t) A(q) = 1 - 0.8109 (+-0.04829) q^-1 C(q) = 1 + 1.231 (+-0.0689) q^-1 + 1.111 (+-0.09903) q^-2 + 0.8187 (+-0.09637) q^-3 + 0.3131 (+-0.06364) q^-4 2
http://www-personal.buseco.monash.edu.au/ hyndman/TSDL/
4
Trayectoria
fac muestral
4
n o i t a l e r r o c o t u A e l p m a S
2 0 −2 −4
0
s n o i t 1 a l e r r o 0.5 c o t u A 0 l a i t r a P −0.5 e l p m −1 a 0 S
100
200
1
0.5
0
−0.5
300
0
10
0
10
facp muestral
20 Lag Variograma
30
30
20
10
10
20
0
30
20
30
Lag
Figura 1: An´alisis de la Serie entrada de gas metano en un horno a gas cu.ft/min El residuo ´o error, indicado en Matlab con e(t) corresponde a un ruido blanco Z n con desviaci´on est´andar estimada dada por σ = 0,1845. La comprobaci´on de que los residuos e(t) son ruido blanco se realiza con la prueba Ljung-Box. En este caso no se rechaza la hip´otesis de incorrelaci´on de e(t). Adem´as, todos los coeficientes resultan significativamente diferentes de cero. Los n´umeros que aparecen entre par´entesis con +- son las desviaciones est´andar sφj de los coeficientes φ j . En este caso, como se explic´o, se puede ver que la dividir el coeficiente por su desviaci´on se obtiene un valor mayor de 1.96 en valor absoluto por lo que se puede considerar que todos son significativamente diferentes de cero.
tcrit = -16.7930 17.8598 11.2191 8.4959 4.9194 Prueba Ljung-Box 0 0
0.1996 0.3949
25.0477 26.2398
31.4104 37.6525
5
Trayectoria
fac muestral
1
n o i t a l e r r o c o t u A e l p m a S
0.5 0 −0.5 −1
0
s n o i t 1 a l e r r o c 0.5 o t u A l a i t r 0 a P e l p m−0.5 a 0 S
100
200
1
0.5
0
−0.5
300
0
10
20 Lag Variograma
facp muestral
30
1.5
1
0.5
10
20
0
30
0
10
20
30
Lag
Figura 2: an´alisis de Residuos del modelo arma(1,4)
2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 −2.5 −3 140
160
180
200
220
240
260
280
Figura 3: an´alisis de pron´osticos con el modelo arma(1,4) 6
300
