Vektorautoregressive Modelle von Produktionsprozessen · J.-F. Emmenegger, D....

Universität Fribourg, Quantitative Wirtschaftsforschung

Vektorautoregressive Modelle vonProduktionsprozessenJean-François Emmenegger1, David Cajander2

1Department für quantitative Wirtschaftsforschung,Universität Freiburg, Schweiz2Holcim Group, [email protected], [email protected]

Schweizer Statistiktage, 25. Oktober 2011, Freiburg, Schweiz

J.-F. Emmenegger, D. Cajander VAR-Simulationsmodelle von HOLCIM-Zementmühlen


Ausgangslage (1)

Die Schweizerfirma HOLCIM betreibt auf 5 KontinentenZementmühlen.Diese Zementmühlen zerkleinern Steingut und Mineralien zukleinkörnigem Sand, mit dem Zement fabriziert wird.Im Frühling 2008 : Herr Ing. D. Cajander, Firma HOLCIM,kontaktiert Herrn Dr. Chr. Leuenberger, Hochschule fürTechnik und Architektur, Freiburg, betreffendZementmühlefragen.Kontakt von Herrn Dr. Chr. Leuenberger mit Dr. J.-F.EmmeneggerEs gibt einen Fragenkomplex um die Zementmühlen.Zementmühlen produzieren Zeitreihen !Im Mai 2008 : El. Ing. M. Bruckhaus und Ing. D. Cajander,Firma HOLCIM, präsentieren umfassend das Problem



Ausgangslage (2)

Erste Vorschläge : Sprünge in den Daten bereinigen, Spektral-und Kointegrations-AnalyseMein Vorschlag : VAR, VARMA, PCA-AnalyseSeptember 2008 : Neue Zusammenkunft in Freiburg :Explizite Vorstellung der Daten und des ProblemsZur Funktionsweise von Zementmühlen kennt man keine DGLoder DiffGL, aber die Maschinen produzieren Zeitreihen !Vorschlag von JFE : Die Zementmühlen aufgrund derZeitdaten durch

1 VAR-Modelle simulieren,2 daraus Eigenschaften der Maschinen ablesen3 und Störverhalten analysieren.



Vektorautoregressive (VAR) Modelle (3)

Man verbindet simultan alle verfügbaren Zeitreihen !Die univariate Methodologie wird erweitert, und zwar vorallem mit dem Ziel ’exogene’ von ’endogenen’ Zeitreihen zuunterscheiden.Liegt zwischen Zeitreihen ’feedback’ vor, so sind vectorautoregressive systems (VAR) zur Beschreibung geeignet.Ist ein VAR-System zu einer Menge von Zeitreihen validiert,so können anschliessend allgemein Impuls-, Interventions-,oder Kausalitätsanalysen angesetzt werden.

’Feedback’ beschreibt die Situation, wenn Informationen übervergangene Resultate Einflüsse auf das gleiche Phänomen in derGegenwart und in der Zukunft haben.



Messdaten von Zementmühlen (4)

0 500 1000 15000

5001000

0 500 1000 15000

500010000

0 500 1000 15000

5000

0 500 1000 15000

5001000

0 500 1000 150050

100150

0 500 1000 15000

100200

0 500 1000 15000

200400

0 500 1000 15000

1020

0 500 1000 15000

100200

0 500 1000 15000

5001000

0 500 1000 15000

5001000

0 500 1000 150050

100150

0 500 1000 15000

1020

0 500 1000 15000

200400

0 500 1000 150050

100150

0 500 1000 15000

50

0 500 1000 15000

50

0 500 1000 15000

50100

0 500 1000 15000

5001000

Figure: 12 ZeitreihenJ.-F. Emmenegger, D. Cajander VAR-Simulationsmodelle von HOLCIM-Zementmühlen



Definitionen Man betrachtet eine Menge diskreterabstandsgleicher Zeitpunkte I , sowie K Zeitreihenwerte xjt ,j = 1, ...,K , die zu einem Vektor (K × 1) Vektorxt = [x1t , x2t , ...., xKt ]′, t ∈ I , zusammengefasst werden.Ferner hat man einen konstanten (K × 1) Vektorνt = [ν1, ν2, ...., νK ]′ und p (K ×K ) Matrizen Ai ,i = 1, ..., p. Ferner hat man ein K -dimensionales WeissesRauschen, auch Innovationsprozess genannt,ut = [x1t , u2t , ...., uKt ]′, mit den Eigenschaften E(ut) = 0,E(utu′t) = Σu und E(utu′s) = 0, für s 6= t. Es wirdvorausgesetzt, dass die Kovarianzmatrix Σu nicht singulär ist,(Luetkepohl).




Damit kann ein vektorautoregressives Modell der Ordnung p,VAR(p) definiert werden.

xt = νt + A1xt−1 + A2xt−2 + ....+ Apxt−p + ut ; t ∈ I (1)

Eine wichtige Bemerkung muss hier angebracht werden. BeiARIMA Modellen ist es notwendig, dass die Zeitreihen stationärsind. Bei der Modellierung von VAR-Systemen ist dies nicht mehrnotwendig, so dass wir vorschlagen können, dass die anfallendenunbereinigten Zeitreihen zu VAR-Systemen modelliert werdenkönnen.



1. Bsp : Simulation eines VAR(1) von 3 Zeitreihen (7)

yt =

543

+ 0.5 0 00.1 0.1 0.3

0 0.2 0.3

yt−1 + u1u2

u3

(2)Erste Simulation mit grossem Rauschen u1 ∼ N (0, 4),u2 ∼ N (0, 4), u3 ∼ N (0, 1).

Simulation of a 3-dim VAR(1) model

VAR(1) with variances:Var(y1)=1,Var(y2)=4),Var(y3)=4

Time Series y1, starting value y1(1)=5

100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y1


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y2


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y3



1. Bsp : Wenig und kein Rauschen im VAR(1) Modell (8)

Zweitens, Simulation mit kleinem Rauschen und ohne Rauschenu1 ∼ N (0, 0.1), u2 ∼ N (0, 0.1), u3 ∼ N (0, 0.1)


VAR(1) with variances:Var(y1)=0.1,Var(y2)=0.1),Var(y3)=0.1


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y1


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y2


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y3


VAR(1) with variances:Var(y1)=0,Var(y2)=0),Var(y3)=0


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y1


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y2


100 200 300 400 500 600 700 800 900 10000.0

2.5

5.0

7.5

10.0

12.5

15.0Y3

y1t = 5 + 0.5y1(t−1)y2t = 4 + 0.1y1(t−1) + 0.1y2(t−1) + 0.3y3(t−1)y3r = 3 + 0.2y2(t−1) + 0.3y3(t−1)

(3)



1. Bsp : Konvergenzen des VAR(1) von 3 Zeitreihen (9)

Mit den Startwerten y11 = 5, y21 = 4 y31 = 3 ergeben sichfolgende Grenzwerte :

y12 = 5 + 0.5y11 = 5 + 0.5× 5 = 7.5y13 = 5 + 0.5y12 = 5 + 0.5× 7.5 = 8.375y1n = 5 + 0.5[1 + 0.5 + 0.52 + ...+ 0.5n ]→ 10

y2t t→∞ = 6.491228

y3t t→∞ = 7.719298

(4)

Dies zeigt deutlich, dass diese Zeitreihen Stabilität oderStationarität aufweisen.Stabilität impliziert Stationarität (Lütkepohl) !



Stabilität von VAR(p) und VAR(1) Modellen (10)

VAR(p)

yt = ν + A1yt−1 + A2yt−2 + ....+ Apyt−p + ut ; t ∈ I (5)

Yt =

[yt

yt−1..

yt−p+1

]; ν =

[ν0..0

]; A =

A1 A2 .. .. Ap−1 ApIK O .. .. O OO IK .. .. O O.. .. .. .. .. ..O O .. .. IK O

; Ut = [ utO..O

].

(6)

det(IKp −Az) = det(Ik −A1z − ...−Apzp) 6= 0 for |z| ≤ 1. (7)

VAR(1)yt = ν + A1yt−1 + ut ; t ∈ I . (8)

det(IK −A1z) 6= 0 for |z| ≤ 1. (9)

The inverse of the roots of the polynomial are the eigenvalues ofthe Matrix !



1. Bsp : Stabilität des VAR(1) Modelles (11)

yt =

543

+ 0.5 0 00.1 0.1 0.3

0 0.2 0.3

yt−1 + u1u2

u3

(10)Stability. The roots of this polynomial are easily seen to bez1 = 2, z2 = 2.1525, z3 = −15.4858. They are obviously all greaterthan 1 in absolute value.Therefore the process (10) is stable. On the other hand the inverseof the roots of the polynomial are the eigenvalues of the Matrix.They have been calculated with RATS, see Doan, and areλ1 = 0.5, λ2 = 0.46458, λ3 = −0.06458. They are smaller than 1.This is the other way to formulate the condition of stability.



1. Bsp : Schätzung des VAR(1) aus Simulation (12)Umkehrung : Aus den 3 simulierten Zeitreihen y1, y2, y3 wird einVAR(p) geschätzt. Das Resultat ist abhängig vom Rauschen !Man kennt die Dimension K = 3.Zuerst muss die Ordnung p des Modells geschätzt werden : AkaikéKriterium von 1973 (AIC) und dem Schwarz BayesianischenKriterium von 1978 (SBC) (RATS, Estima, Evanston, Ill, USA.).AIC und SBC indizieren beide : lag p = 1Elimination der nicht sign. Matrixkoeffizienten -> Near-VAR !

i j var. DW σ̂ aij const. Std. err. p-value(1) (2) (3) (4) (5) (6) (7) (8) (9)1 1 y1 2.0381 0.1011 0.51749 0.0242 0.0000

4.8151 0.2423 0.00002 1 y2 1.9872 0.1010 0.09675 0.0251 0.00012 2 0.12000 0.0290 0.00002 3 0.30998 0.0292 0.0000

3.8181 0.3076 0.00003 2 y3 2.0283 0.0961 0.18627 0.0277 0.00003 3 0.31997 0.0283 0.0000

2.9753 0.2402 0.0000

Table: Schätzung des Near-VAR(1) Modelles (ohne Rauschen)



1. Bsp : Diagnose des gefundenen VAR(1) Modelles (13)

Die geschätzten Koeffizienten werden mit den ursprünglichenKoeffizienten verglichen.Dies geschieht separat für jede Gleichung des VAR(1) Modelles

yt = ν + A1yt−1 + ut ; t ∈ I . (11)

i j var. const. ĉonst. âij aij F-test p-value(1) (2) (3) (4) (5) (6) (7) (8) (9)

y1 5 4.8151 0.5819 0.44571 1 0.51749 0.5 0.5184 0.4717

y2 4 3.81897 0.3497 0.55442 1 0.09675 0.1 0.0166 0.89732 2 0.12000 0.1 0.0290 0.49002 3 0.30998 0.3 0.0292 0.7327

y3 3 2.9753 0.2402 0.91833 2 0.18627 0.2 0.0277 0.62043 3 0.31997 0.3 0.0283 0.4801

Table: Joint F-Wald tests on the VAR(1)-model coefficients

Je grösser das Rauschen, desto schlechter die Uebereinstimmung !



1. Bsp : Impulsantworten des VAR(1) Model (14)

s1 = [1, 0, 0]′; s2 = [0, 1, 0]′; s3 = [0, 0, 1]′ (12)eingesetzt in die Startwerte : y0 = si , i=1,2,3

yt =

0.5 0 00.1 0.1 0.30 0.2 0.3

yt−1 + ut ; t = 1, 2, 3, ... (13)mit verschiedenen Fehlervektoren : ut = [u1t , u2t , u3t ]′. Es gibtkeinen Konstantenvektor, da ein Impuls ein additives Signal ist !

t s1 s2 s3(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)0 1 0 0 0 1 0 0 0 11 0.5 0.1 0 0 0.1 0.2 0 0.3 0.32 0.25 0.06 0.02 0 0.07 0.08 0 0.12 0.153 0.125 0.037 0.018 0 0.07 0.038 0 0.057 0.069

Table: Impulsantworten zu s1, s2, s3, (ohne Rauschen)




Responses of each series yi to each of the initial shocks

IR of y1 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (1,1)

impblk (1,2)impkbl (1,3)

IR of y2 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (2,1)


IR of y3 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (3,1)


Figure: Zeitreihen Antworten auf Impulse : VAR(1), (1)




Intervention Analysis with a VAR(1) Model

Impulse Responses (IR) to initial shocks at time t=0

IR of y1, y2, y3 to shock [1,0,0]

5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (1,1)impblk (2,1)impkbl (3,1)


5 10 15 200.00

0.25

0.50

0.75

1.00

1.25


IR of y1, y2 ,y3 to shock [0,0,1]

5 10 15 200.00

0.25

0.50

0.75

1.00

1.25


Figure: Impuls Antworten der Zeitreihen : VAR(1), (1)




Responses of each series yi to each of the initial shocks

IR of y1 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (1,1)


IR of y2 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (2,1)


IR of y3 to shocks [1,0,0], [0,1,0], [0,0,1]

0 5 10 15 200.00

0.25

0.50

0.75

1.00

1.25

1.50impblk (3,1)


Figure: Zeitreihen Antworten auf Impulse : VAR(1), (2)




Intervention Analysis with a Near-VAR(1) Model

Fig. 3 Impulse Responses (IR) to initial shocks at time t=0


5 10 15 200.00

0.25

0.50

0.75

1.00

1.25



5 10 15 200.00

0.25

0.50

0.75

1.00

1.25


IR of y1, y2 ,y3 to shock [0,0,1]

5 10 15 200.00

0.25

0.50

0.75

1.00

1.25


Figure: Impuls Antworten der Zeitreihen : VAR(1), (2)



2. Bsp : VAR Modell für A2 und P2 (19)

xt = νt +24∑

i=1Aixt−i + ut ; t ∈ I . (14)

Es wurden in eimem ersten VAR-Modellversuch mit K = 2 diebereingten Zeitreihen Speed main drive A2t und Torque peakto Peak P2t ausgewählt, t ∈ I .Das erste Ziel der VAR-Modellierung besteht darin, zumgegebenen Zeitreihenvektor xt = [A2t ,P2t ]′ den GradMaxlag=p zu bestimmen und entsprechend das Modell mitmöglichst wenig Parametern zu determinieren (parsimoniousmodel).




Bereinigte Holcim-Series A2, P2

100 200 300 400 500 600 700 800 900 1000 1100 1200-3.0

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

A2

100 200 300 400 500 600 700 800 900 1000 1100 1200-8

-6

-4

-2

0

2

4

6

8

P2

Figure: 2-dim Zeitreihen-Modell




In diesem Sinne wird zuerst ein VAR(p)-Modell vomvorläufigen Grade p = 24 angesetzt.

xt = νt +24∑

i=1Aixt−i + ut ; t ∈ I . (15)

Ein solches sog. strukturiertes VAR(p)-Modell ist natürlichüberparametrisiert ; es hat nämlich n = K + p ·K 2 = 98Parameter.Bestimmung des optimalen Grades mit dem AkaikeInformation Criterion (AIC) und dem Schwarz-Bayes Criterion(SBC).Das AIC : p = 13, SBC : p = 8. Man wählt den Grad p = 13.




Elimination aller nicht signifikanten Matrixelemente derMatrizen Ai , i = 1, .., 13.Der Konstantenvektor ν ist nicht signifikant.Es verbleiben noch n = 14 Parameter im nun entstandenennear-VAR Modell.

A1 =[

0.162 0.9490.141 0

]; A2 =

[0 0.0800 0

]; A3 =

[0 0

0.123 0

]A4 =

[0 00 0

]; A5 =

[0 00 0

]; A6 =

[0 −0.529

0.160 0

]A7 =

[0 0.4970 0

]; A8 =

[0.100 00.100 0

]; A9 =

[0 00 0

]A10 =

[0 00 −0.016

]; A11 =

[0 00 0.025

]; A12 =

[0 −0.2730 0

]A13 =

[0 00 0.245

]Wenige Parameter im VAR-Modell, keine Gefahr derUeberparametrisierung.




Ein near-VAR Modell wird nun mit den obigen 13 Matrizengebildet :

xt =13∑

i=1Aixt−i + ut ; t ∈ I . (16)

Es müssten noch die Residuen untersucht werden.Dann wird der Likelihood ratio test von Sims angewandt, umdie Verbesserung der Modellierung beim Übergang vomstrukturierten VAR(24) (15) zum parsimoniertenVAR(13)-Modell (16). Der χ2-Test spricht an, und dieausgewiesene Statistik ist P(χ2(40) ≥ 36.02) = 0.6500.Die Verbesserung ist somit offensichtlich.



2. Bsp : Companion-Matrix A (24)

A =

0.16 0.95 0 0.08 0 0 ... 0 0 0 0 0 0 −0.27 0 00.14 0 0 0 0.12 0 ... 0 0 −0.016 0 0.03 0 0 0 0.25

1 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 1 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 1 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 1 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 1 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 1 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 1 0 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 1 0 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 1 0 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 1 0 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 1 0 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 1 0 0 00 0 0 0 0 0 ... 0 0 0 0 0 0 1 0 0

(17)



2. Bsp : Eigenwerte der Matrix A (25)

A =

A1 A2 .. .. A12 A13I2 O .. .. O OO I2 .. .. O O.. .. .. .. .. ..O O .. .. I2 O

(18)

i e(i) i e(i) i e(i) i e(i) i e(I)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)1 0.95799 7 -0.69162 13 -0.55002 19 0.09616 25 -0.000322 -0.92684 8 -0.69162 14 -0.53862 20 0.01036 26 0.000203 -0.92684 9 -0.69162 15 0.43686 21 0.010364 0.85899 10 0.58946 16 -0.31079 22 0.000855 0.78825 11 0.58946 17 -0.31079 23 -0.000726 0.78825 12 -0.55002 18 0.09616 24 -0.00065

Table: Eigenvalues of matrix A (17)



3. Bsp : VAR Modell A1, P1, Q1 (26)

Raw Holcim-Series A1, P1, Q1A1

100 200 300 400 500 600 700 800 900 1000 1100 1200830.0

832.5

835.0

837.5

840.0

842.5

845.0

847.5

P1

100 200 300 400 500 600 700 800 900 1000 1100 12000

10

20

30

40

50

Q1

100 200 300 400 500 600 700 800 900 1000 1100 12000

10

20

30

40




3. Bsp : VAR Modell für A1, P1, Q1 (27)

Ein VAR-Modell mit den K = 3 nicht bereinigten ZeitreihenA1t , P1t und Q1tBildung des Zeitreihenvektors xt = [A1t ,P1t ,Q1t ]′ vom GradpBestimmung des Grades p mit AIC und SBC KriteriumEine RATS Routine berechnet p = 13 mittels AIC und p = 7mit SBC.Ansatz mit VAR(20) gegen VAR(13) : Likelihood ratio testvon Sims P(χ2(63) ≥ 37.295) = 0.9959, also VAR(13)p = 13 kann nicht weiter herunter gedrückt werden !Hingegen : Ein Near −VAR(13) Modell mit nurmehr total 20verschiedenen Parametern ist validiert, wobei der Likelihoodratio test von Sims das reduzierte Modell Near −VAR(13)nur knapp annimmt : P(χ2(99) ≥ 122.531) = 0.05456.



3. Bsp : VAR Modell für A1, P1, Q1 (28)

xt = νt +13∑

i=1Aixt−i + ut ; ut = i.i.d.(0, Σ̂u); (19)

Autocorrelations of residuals of a VAR(13) model

Fig. x VAR(13) model with series A1, P1, Q1

A1 and A1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and P1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and Q1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Autokorrelation : Keine serielle Korrelation mehr !Portemanteau Test von Hosking auf VAR(10) :P(χ2(153) ≥ 136.2133) = 0.8311.



3. Bsp : VAR Modell A1, P1, Q1 und ImpRA (29)

Bisherige Ergebnisse zu den nichtbereinigten Zeitreihen A1t ,P1t und Q1t heranziehen.Man nimmt also das VAR(10) und das VAR(13) Modell.Holcim : die Zeitreihe A1t ist signalverursachend (exogen)die Zeitreihen P1t und Q1t wiedergeben den Einfluss dersignalverursachenden (endogen) Zeitreihen.Ansatz zum Zeitpunkt t = 0 :X1 = [A1(0) = 1,P1(0) = 0,Q1(0) = 0]′,X2 = [A1(0) = 0,P1(0) = 1,Q1(0) = 0]′,X3 = [A1(0) = 0,P1(0) = 0,Q1(0) = 1]′.Ausgehend von diesen drei Impulsvektoren X1, X2, X3 inVAR(13) und VAR(10) Impulsantworten zu zukünftigenZeitpunkten t > 0 analysieren.



3. Bsp : VAR Modell A1, P1, Q1 und ImpRA (30)

Impulsvektoren zum Zeitpunkt : t = 0X1 = [A1(0) = 1,P1(0) = 0,Q1(0) = 0]′,X2 = [A1(0) = 0,P1(0) = 1,Q1(0) = 0]′,X3 = [A1(0) = 0,P1(0) = 0,Q1(0) = 1]′.

Intervention Analysis with VAR(13)


IR to shock [1,0,0]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

IR to shock [0,1,0]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

IR to shock [0,0,1]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

Intervention Analysis with VAR(10)


IR to shock [1,0,0]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

IR to shock [0,1,0]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

IR to shock [0,0,1]

5 10 15 20-1.0

-0.5

0.0

0.5

1.0

1.5A1

P1

Q1

Figure: Impulsantworten



4. Bsp : VAR für A1, P1, Q1, B1, C1, D1 (31)

Rohe Holcim-Series A1, P1, Q1, B1, C1 D1

250 500 750 1000830.0

832.5

835.0

837.5

840.0

842.5

845.0

847.5

A1

250 500 750 10000

10

20

30

40

50

P1

250 500 750 10000

5

10

15

20

25

30

35

40

45

Q1

250 500 750 10000

1000

2000

3000

4000

5000

6000

7000

B1

250 500 750 10000

500

1000

1500

2000

2500

3000

3500

4000

4500

C1

250 500 750 1000200

300

400

500

600

700

800

D1




4. Bsp : VAR für A1, P1, Q1, B1, C1, D1 ; ImpRA (32)

Analyse der Zeitreihen A1, P1, Q1, B1t , C1t , D1tAIC und BSC Informationskriterium : Verzögerung p = 8Man bestimmt also ein strukturiertes VAR(8).Portemanteau test von Hosking :P(χ2(504) ≥ 458.8265) = 0.9260. Okay !Holcim : die Zeitreihe A1t ist Signal verursachend, dieZeitreihen P1t , Q1t , B1t , C1t , D1t geben den Einfluss derSignal verursachenden Zeitreihen wieder.Impulsvektoren zum Zeitpunkt t = 0 :X1 = [A1(0) = 1,P1(0) = 0,Q1(0) = 0,B1(0) = 1,C1(0) = 0,D1(0) = 0]′X2 = [A1(0) = 0,P1(0) = 1,Q1(0) = 0,B1(0) = 0,C1(0) = 0,D1(0) = 0]′X3 = [A1(0) = 0,P1(0) = 0,Q1(0) = 1,B1(0) = 0,C1(0) = 0,D1(0) = 0]′X4 = [A1(0) = 0,P1(0) = 0,Q1(0) = 0,B1(0) = 1,C1(0) = 0,D1(0) = 0]′X5 = [A1(0) = 0,P1(0) = 0,Q1(0) = 0,B1(0) = 0,C1(0) = 1,D1(0) = 0]′

X6 = [A1(0) = 0,P1(0) = 0,Q1(0) = 0,B1(0) = 0,C1(0) = 0,D1(0) = 1]′



4. Bsp : VAR für A1, P1, Q1, B1, C1, D1 ; ImpRA (33)

Autocorrelations of residuals of a VAR(8) modelA1 and A1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and B1

0 5 10 15

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and C1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

A1 and D1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and P1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and B1

0 5 10 15

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and C1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

P1 and D1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and Q1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and B1

0 5 10 15

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and C1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Q1 and D1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and B1

0 5 10 15 20

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and C1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

B1 and D1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and B1

0 5 10 15

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and C1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

C1 and D1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and A1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and P1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and Q1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and B1

0 5 10 15

-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and C1

0 5 10 15-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

D1 and D1

0 5 10 15 20-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00

Figure: Autokorrelationen der VAR(8)-Vektorresiduen



4. Bsp : VAR für A1, P1, Q1, B1, C1, D1 ; ImpRA (34)Intervention Analysis

IR to shock [1,0,0,0,0,0]

5 10 15 20 25 30-75

-50

-25

0

25

50

75A1

P1

Q1

B1

C1

D1

IR to shock [0,1,0,0,0,0]

5 10 15 20 25 30-40

-30

-20

-10

0

10

20

30

40A1

P1

Q1

B1

C1

D1

IR to shock [0,0,1,0,0,0]

5 10 15 20 25 30-40

-30

-20

-10

0

10

20

30

40A1

P1

Q1

B1

C1

D1

Figure: Impulsantworten von A1, P1, Q1 auf Anfangsschocks



4. Bsp : VAR für A1, P1, Q1, B1, C1, D1 ; ImpRA (35)Intervention Analysis

IR to shock [0,0,0,1,0,0]

5 10 15 20 25 30-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00A1

P1

Q1

B1

C1

D1

IR to shock [0,0,0,0,1,0]

5 10 15 20 25 30-5.0

-2.5

0.0

2.5

5.0A1

P1

Q1

B1

C1

D1

IR to shock [0,0,0,0,0,1]

5 10 15 20 25 30-20

-15

-10

-5

0

5

10

15

20A1

P1

Q1

B1

C1

D1

Figure: Impulsantworten von B1, C1, D1 auf Anfangsschocks



4. Bsp : Impulsantwort Analyse (36)Intervention Analysis

IR to shock [1,0,0,0,0,0]

5 10 15 20 25 30-1.00

-0.75

-0.50

-0.25

0.00

0.25

0.50

0.75

1.00A1P1

Q1

Figure: Impulsantworten von B1, C1, D1 auf Anfangsschock[1,0,0,0,0,0]’



5. Bsp : Signale der Dataset21 (37)

Exogene SignaleSMDpath : Speed Main Drive ist die Motor-Geschwindigkeit der Mühle (Mahlteller)PHpath : Pressure Hydraulic ist der Druck (in bar) der Mahlroller in der Mühle.SpSepPath : Speed Separator ist die Geschwindigkeit des Separators.FPpath : Flow Pump ist die Wasser-Injektion. Dies hilft das Material Mahlbett kompakter zu machen, sodass die Mühle stabiler läuft.TT3path : Temperatur T3 ist die Messung der Gastemperatur am Mühleausgang (TT3 könnte endogenund Temp12Average exogen sein).GFpath : Gas Flow ist der Gas-Fluss durch die Mühle und ist durch einen grossen Ventilator kontrolliert.VRMSRpath : Vertical Roller Mill Slag Rate ist die Materialeingabe in die Mühle

Endogene SignaleCurrentMainDrive : Strom im Hauptmotor.PositionRoller : Messung der Position der Mahlroller, dies gibt eine Information über die Materialbetthöheauf dem Mahlteller.Temp12Average : Eingang der Gastemperatur. Da es zwei Eingangsmessungen T1 und T2 gibt, hat maneinen Durchschnittswert gerechnet.Vibration : Vibration der MühleTorquePtP : Drehmoment Peak To Peak auf Motorachse gemessen.



5 Bsp : Eingangs-Dataset21 (38)

AsymVAR Forecasts (VRM Short Model/Scaled Steps/Dataset21)Speed Main Drive [rpm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500839.5

840.0

840.5

841.0

841.5

842.0

842.5

843.0

Pressure Hydraulic1 [bar]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50066

67

68

69

70

71

GasFlow [m3/h]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500390000

410000

430000

450000

470000

TempT3 [degC]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50080

85

90

95

100

105

Speed Separator [rpm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500120.0

125.0

130.0

135.0

140.0

Flow Pump [l/h]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500700

800

900

1000

1100

1200

1300

VRM slag rate [t/h]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50080

90

100

110

120

130

140

Current Main Drive [A]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500600

650

700

750

800

Position Roller [mm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50025

35

45

55

65

Vibration [mm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50011

13

15

17

19

Temp 1&2 Average [degC]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500160

170

180

190

200

TorquePtP [kNm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50022

24

26

28

30

32

34

36




RATS-Analyse von Dataset21 (39)

Visualisierung von Dataset21Visualisierung der Pfad-Zeitreihen und endogenen Zeitreihen(*)Definition of Asymmetric (Near) VAR with MaxLag=8, usingthe RATS-Routine ’stwise’ to eliminate all non-significantmatrix elementsIntroduction of the Intervention-Analyses with the pathesCalculate the in-sample-forecasts on the endogenes series onthe basis of (*)



5. Bsp : Struktogram des Algorithmus (40)

i e(i) i e(i) i e(i) i e(i) i e(I)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)1 0.95799 7 -0.69162 13 -0.55002 19 0.09616 25 -0.000322 -0.92684 8 -0.69162 14 -0.53862 20 0.01036 26 0.000203 -0.92684 9 -0.69162 15 0.43686 21 0.010364 0.85899 10 0.58946 16 -0.31079 22 0.000855 0.78825 11 0.58946 17 -0.31079 23 -0.000726 0.78825 12 -0.55002 18 0.09616 24 -0.00065

Tabelle 5: Eigenvalues of matrix A (40)

Here exogenous (=Input) means that the corrresponding signals are influencedby the driving operation of the machine, whereas endogenous (=output) means thatthe signal is generated by the machine and is the result of the action on the exogenoussignals.

For each of the 5 endogenous time series a linear regression meodel is created,depending on all the 12 variables. The maximal lag is set to lag = 8. The process isa so-called stepwise regression, taking only the significant lagged coefficients of eachof the the depending variables.

Then a Near-VAR model is created with these 5 linear regression time seriesand the Near-VAR model gets the name: Nearvar.

As last step, the in sample forecasts of the 7 exogenous variables are calculatedover a horizon of 400 values. STRUKTOGRAM OF THE MAIN NEAR-VAR PROGRAM

******************************************************open data files with 12 initial series

SpeedMainDrive CurrentMainDrive PressureHydraulic1 PositionRoller GasFlow TempT3TorquePtP SpeedSeparator FlowPump Vibration Temp12Average VRMslagrate

open data files with 5 initial endogenous PATH series and 7 exogenous seriesSMDpath CurrentMainDrive PHpath PositionRoller GFpath TT3path TorquePtPSpSepPath FPpath Vibration Temp12Average VRMSRpath

* NEAR VAR Model of 5 endogenous time series with constant and all the 12 time series*STEP WISE procedure, with auto tuning, with lag = 8

stwise(define=CMDEq) CurrentMainDrivestwise(define=PREq) PositionRollerstwise(define=VibEq) Vibrationstwise(define=T12AEq) temp12Averagestwise(define=TPtPEq) TorquePtP

*Create a System and estimatesystem CMDEq PREq VibEq T12AEq TPtPEqestimate(print,residuals=resvar)

*Group the endogenous tiem seriesGroup NearVAR CMDEq PREq VibEq T12AEq TPtPEq

*The exogenous variables are replaced by the path values(for the step):set SpeedMainDrive = SMDpathset PressureHydraulic1 = PHpathset GasFlow = GFpathset TempT3 = TT3pathset SpeedSeparator = SpSepPathset FlowPump = FPpathset VRMslagrate = VRMSRpath

*Forecastforecast(model=NearVAR,results=forecasts,from=100,to=502)

*Graph and print times series with forecastSPGRAPH the 12 time seriesPrint forecasts

* Create an EXCEL Dataset

Datasetforecast.xls

21



5. Bsp : In-sample forecast mit Dataset21 (41)

T-Stat/Sig Auto-tuned Near/AsymVAR Forecasts (VRMSR,DS21)SMDpath [rpm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500839.5

840.0

840.5

841.0

841.5

842.0

842.5

843.0

PHpath [bar]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50066

67

68

69

70

71

GFpath [m3/h]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500390000

410000

430000

450000

470000

TT3path [degC]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50080

85

90

95

100

105

SpSepPath [rpm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500120.0

125.0

130.0

135.0

140.0

FPpath [l/h]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500700

800

900

1000

1100

1200

1300

VRMSRpath [tph]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50060

70

80

90

100

110

120

130

FORECAST Current Main Drive [A]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500600

650

700

750

800

FORECAST Position Roller [mm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50025

35

45

55

65

FORECAST Vibration [mm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50011

13

15

17

19

FORECAST Temp12Average [degC]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 500150

160

170

180

190

200

210

FORECAST TorquePtP [kNm]

25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425 450 475 50022

24

26

28

30

32

34

36




Gegenwärtiger Stand (42)

Algorithmus unabhängig von einer lizenzierten Softwareumschreiben !Online Applikation aufbauen, testen und einführen !Quantitative Verwendung !

14

© Holcim Group Support Limited 2009 Remote Diagnostic – 20.04.2009/CAJConfidential – reserved to Process Innovation Committee and PIP Team

Details Graphical Interpretation: Flow Apron Feeder

The feed increase is clearly reflected on the current, power and torque signals.

All pth lag model forecasts are consistent: this demonstrates that the equations parameters could be correctly estimated.

Quantitative verification of ∆P(∆I): PActive= Ieff·Ueff·cos(φ)·31/2

- with cos(φ)=0.85 and Ueff=4kV

- PCalc,t1=4’300kW → PFcast,t1= 4’190kW- PCalc,t2=4’390kW → PFcast,t2= 4’280kW

RelativeError(PActive) = ∆P/P = 2.5%

Given the incertitude on cos(φ) and Ueff the relative error is acceptable and demonstrate strong quantitative accuracy of the forecasts.

The proportionality between current and power is conserved.

Analysis Results7.



Vektorautoregressive Modelle von Produktionsprozessen · J.-F. Emmenegger, D....

Documents

Transcript of Vektorautoregressive Modelle von Produktionsprozessen · J.-F. Emmenegger, D....