Procesos Estocásticos

Taller No 2: Análisis procesos ARMA con Matlab

Octubre 2007

1. Funciones Matlab para análisis de procesos AR-

MA
1. Funciones para identificación de los órdenes p y q:

a) autocorr: calcula y grafica la función de autocorrelación estimada.

b) parrcorr: calcula y grafica la función de autocorrelación parcial estimada.

2. Funciones para estimación y ajuste:

a) armabat: función para identificar la pareja (p,q) que produce el modelo

con menor criterio de información de Akaike (AIC). Ver pag. web de H.
Hurd. 1
b) armax: función que estima los parémetros del modelo ARMA(p,q). Mat-
lab provee varias funciones para estimación de modelos pero solamente se
considerará ésta.
c) resid: función para calcular los residuos del modelo con el fin de poder
realizar pruebas de hipótesis para determinar si es ruido blanco, como la
prueba de Ljung-Box. Los residuos son estimaciones de Zn .
d ) lbqtest: función para realizar la prueba de Ljung-Box.
e) compare: función para examinar la calidad de los pronósticos que se pueden
hacer con el modelo ajustado con el fin de determinar si éste es adecuado
o nó.

Notación Matlab. En lugar de la letra L, Matlab utiliza q −1 , luego q −1 (Xn ) = Xn−1 .

Por ejemplo, un modelo ARMA(4,2) se expresa en Matlab ası́:
A(q −1)Xn = B(q −1)Zn
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álisis con las funciones de Matlab

1. Graficar la trayectoria del proceso, la FAC estimada, la FACP estimada y el

Variograma. Las instrucciones a continuación muestran cómo 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és de una posible identificación del tipo de proceso se procede a especificar

los órdenes p y q del proceso. En caso de no ser posible identificar un AR ó un
MA, se toma inicialmente el rango p, q = 1, 2, 3, 4, 5, 6, y se corre la función
“armabat”. Antes de aplicar la función 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án los valores de los órdenes p y q que
mejor describen el proceso.

3. Estimación los parámetros del modelo ARMA(p,q). Para esto se utiliza la fun-
ción “armax” con la pareja (p, q) escogida en el punto anterior. Esta instrucción
crea un objeto de nombre (por ejemplo)“armapq” que contiene varios campos
con información sobre el modelo estimado.

armapq = armax(xt,[pbest qbest]);

present(armapq)
armapq.a
armapq.c

4. Pruebas de significación de los Parámetros. La instrucción

tcrit = armapq.ParameterVector./sqrt(diag(armapq.CovarianceMatrix))

calcula un vector de cocientes tj = φbj /Sφbj , donde Sφbj es la desviación estándar

bj . Cada valor tj es el valor de un estadı́stico t de Student que
del coeficiente φ
sirve para probar la hipótesis de que los coeficientes son significativamente difer-
entes de cero, es decir, se acepta la hipótesis Ho : φj = 0 al nivel de 5 % si
observamos que |tj | < 1,96; en caso contrario se rechaza.
También se puede calcular el valor estimado de la varianza del ruido blanco, σ 2,
mediante la instrucción: armapq.NoiseVariance.

5. Para completar el análisis es necesario chequear si los residuos del modelo ajus-
tado son ruido blanco. Los residuos son valores estimados del proceso Zn . 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álculo 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ó bien los datos es ajustar
el modelo con la primera mitad de los datos y utilizar la parte restante para
comparar con los pronósticos a un paso: se compara Xn con el pronóstico de Xn
realizado con el modelo. La función “compare” de Matlab hace este cálculo.

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áfica 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 par-
cial estimadas. Según lo explicado acerca de la identificación de modelos tipo AR-
MA, corresponderı́a a un proceso autorregresivo AR(3). La estimación con armabat
y armax se realizó de acuerdo a las indicaciones anteriores y se obtuvo un modelo
ARMA(1,4)(!!!). El modelo ajustado para la serie Xn − µ 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 1

Sample Autocorrelation
2
0.5
0
0
−2

−4 −0.5
0 100 200 300 0 10 20 30
Lag
Sample Partial Autocorrelations

facp muestral Variograma

1 30

0.5
20
0
10
−0.5

−1 0
0 10 20 30 0 10 20 30
Lag

Figura 1: Análisis de la Serie entrada de gas metano en un horno a gas cu.ft/min

El residuo ó error, indicado en Matlab con e(t) corresponde a un ruido blanco Zn

con desviación estándar estimada dada por σ b = 0,1845. La comprobación de que los
residuos e(t) son ruido blanco se realiza con la prueba Ljung-Box. En este caso no se
rechaza la hipótesis de incorrelación de e(t). Además, todos los coeficientes resultan
significativamente diferentes de cero. Los números que aparecen entre paréntesis con
+- son las desviaciones estándar sφbj de los coeficientes φ bj . En este caso, como se
explicó, se puede ver que la dividir el coeficiente por su desviación 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.1996 25.0477 31.4104

0 0.3949 26.2398 37.6525

5
 Trayectoria fac muestral
1 1

Sample Autocorrelation
0.5
0.5
0
0
−0.5

−1 −0.5
0 100 200 300 0 10 20 30
Lag
Sample Partial Autocorrelations

facp muestral Variograma

1 1.5

0.5 1

0 0.5

−0.5 0
0 10 20 30 0 10 20 30
Lag

Figura 2: análisis de Residuos del modelo arma(1,4)

1.5

0.5

−0.5

−1

−1.5

−2

−2.5

−3
140 160 180 200 220 240 260 280 300

Figura 3: análisis de pronósticos con el modelo arma(1,4)