1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
|
C Last change: SRD 31 Jan 100 7:39 am
SUBROUTINE setmdl(Estprm,Laumts)
IMPLICIT NONE
c-----------------------------------------------------------------------
c setmdl.f, Release 1, Subroutine Version 1.3, Modified 01 May 1995.
c-----------------------------------------------------------------------
c Subroutine to calculate exact MA ARIMA filter residuals.
c Setmdl differences the X:y matrix and changes the model to remove
c the differencing. The the remaining ARMA model is estimated.
c Model information is in ARIMA.cmn common so the variables are saved
c between calls of the routines fcnkf, and arflt. Setmdl also
c constructs a vector of parameters to be estimated in the nonlinear
c routine. Fcnar calculates ARIMA filter residuals given new estimated
c parameters, estprm, from the nonlinear routine, regression residuals,
c tsrs, from rgcpnt, and the model information that was constructed
c in setmdl. ARflt filters an extended [X:y] matrix from rgcpnt
c using parameter estimates saved during the last fcnkf call.
c-----------------------------------------------------------------------
c Jan 2000 - Argument added to facilitate testing if initial values
c generated from HR routine were valid (BCM)
c-----------------------------------------------------------------------
c Name Type Description
c-----------------------------------------------------------------------
c begptr i Local pointer to the first row in opr of the current
c difference, AR, or MA filter
c estprm d Input nestpm long vector of estimated parameters from the
c nonlinear routine.
c iflt i Local index for the current filter type, DIFF, AR, or MA.
c ilag i Local index for the current lag, pointer to the current
c element in lag,coef, and fix.
c iopr i Local index for the current operator, it is the pointer
c to the current row in the operator specfication matrix,
c opr.
c beglag i Local pointer to the current coefficient and lag in arimap
c and arimal
c nestpm i Input number of parmeters in estprm
c nlag i Local number of lags in the current operator of a filter
c nopr i Local for the number of operators in a DIFF, AR, or MA
c filter.
c one d Local PARAMETER for a double precision 1
c oprptr i Local pointer of the current row of opr, specifying the
c current operator
c zero d Local PARAMETER for double precision 0
c-----------------------------------------------------------------------
INCLUDE 'stdio.i'
INCLUDE 'srslen.prm'
INCLUDE 'model.prm'
INCLUDE 'model.cmn'
INCLUDE 'mdldat.cmn'
INCLUDE 'units.cmn'
INCLUDE 'error.cmn'
c ------------------------------------------------------------------
LOGICAL T,F
PARAMETER(T=.true.,F=.false.)
c ------------------------------------------------------------------
INTEGER beglag,begopr,endlag,endopr,iflt,ilag,iopr
LOGICAL Laumts
DOUBLE PRECISION Estprm
DIMENSION Estprm(PARIMA)
c-----------------------------------------------------------------------
LOGICAL dpeq
EXTERNAL dpeq
c-----------------------------------------------------------------------
c The following declaractions were added by Bor-Chung Chen for
c the initial checking of invertibility and stationarity on 4/26/1995.
c-----------------------------------------------------------------------
CHARACTER dotln*(POPRCR+1),tmpttl*(POPRCR)
LOGICAL allinv,allfix,onunit,mains,maon,mainsa,maonua,arona
INTEGER factor,degree,ntmpcr,i
DOUBLE PRECISION coef(PORDER+1),zeror(PORDER),zeroi(PORDER),
& zerom(PORDER),zerof(PORDER)
c-----------------------------------------------------------------------
c The following declaractions were added by Brian C. Monsell for
c the shrinkage of MA operators when roots are close to the unit circle
C on 1/8/1998.
c-----------------------------------------------------------------------
DOUBLE PRECISION PT9
PARAMETER(PT9=0.9D0)
LOGICAL shrnkp
INTEGER lagind
DIMENSION lagind(PORDER)
LOGICAL first
SAVE first
DATA first/.true./
c-----------------------------------------------------------------------
DATA dotln/
& ' -----------------------------------------------------------'
& /
c-----------------------------------------------------------------------
mainsa=F
maonua=F
arona=F
c-----------------------------------------------------------------------
c For each AR and MA filter find out how many estimated parameters
c and add them to the estimated parameter vector estprm.
c-----------------------------------------------------------------------
Nestpm=0
shrnkp=F
DO iflt=DIFF,MA
begopr=Mdl(iflt-1)
endopr=Mdl(iflt)-1
c ------------------------------------------------------------------
DO iopr=begopr,endopr
beglag=Opr(iopr-1)
endlag=Opr(iopr)-1
c ------------------------------------------------------------------
DO ilag=beglag,endlag
IF(.not.Arimaf(ilag))THEN
Nestpm=Nestpm+1
Estprm(Nestpm)=Arimap(ilag)
c-----------------------------------------------------------------------
c CODE ADDED BY Brian Monsell Jan. 1998
c Save lag corresponding to the ith estimated operator, in case
c this operator is to be shrunk later.
c-----------------------------------------------------------------------
IF(.not.first)lagind(ilag)=Nestpm
c-----------------------------------------------------------------------
END IF
END DO
END DO
END DO
c-----------------------------------------------------------------------
c Compute the roots of initial theta(B)=0
c-----------------------------------------------------------------------
begopr=Mdl(MA-1)
beglag=Opr(begopr-1)
endopr=Mdl(MA)-1
c ------------------------------------------------------------------
IF(endopr.gt.0)THEN
endlag=Opr(endopr)-1
mainsa=F
maonua=F
c ------------------------------------------------------------------
DO iopr=begopr,endopr
beglag=Opr(iopr-1)
endlag=Opr(iopr)-1
c ------------------------------------------------------------------
factor=Oprfac(iopr)
degree=Arimal(endlag)/factor
coef(1)=-1.0D0
CALL setdp(0D0,degree,coef(2))
c ------------------------------------------------------------------
CALL setdp(0D0,PORDER,zeror)
CALL setdp(0D0,PORDER,zeroi)
CALL setdp(0D0,PORDER,zerom)
c ------------------------------------------------------------------
DO ilag=beglag,endlag
coef(Arimal(ilag)/factor+1)=Arimap(ilag)
END DO
c ------------------------------------------------------------------
CALL roots(coef,degree,allinv,zeror,zeroi,zerom,zerof)
c-----------------------------------------------------------------------
c Check invertibility
c the roots are g(i)=(zeror(i), zeroi(i)), i=1,2,...,degree
c complex roots are g(i) and g(i+1)
c If all zeros are invertible do nothing; otherwise print an error
c message and STOP execution of the program.
c-----------------------------------------------------------------------
mains=F
IF(.not.allinv)THEN
c-----------------------------------------------------------------------
CALL getstr(Oprttl,Oprptr,Noprtl,iopr,tmpttl,ntmpcr)
IF(Lfatal)RETURN
IF(.not.Laumts)WRITE(STDERR,1010)tmpttl(1:ntmpcr)
WRITE(Mt2,1010)tmpttl(1:ntmpcr)
1010 FORMAT(/,' ERROR: ',a,' polynomial with initial parameters',
& ' is noninvertible',/,' with root(s) inside the',
& ' unit circle. RESPECIFY model with',/,
& ' different initial parameters.',/)
c ------------------------------------------------------------------
mains=T
mainsa=T
c ------------------------------------------------------------------
END IF
onunit=F
c-----------------------------------------------------------------------
c CODE ADDED BY Brian Monsell Jan. 1998
c Initialize shrnkp, which indicates whether the operator should
c be shrunk.
c-----------------------------------------------------------------------
shrnkp=F
i=0
DO WHILE (i.lt.degree)
i=i+1
c-----------------------------------------------------------------------
c CODE ADDED BY Brian Monsell Jan. 1998
c If first entry, to see if initial parameters are noninvertible
c with roots on the unit circle. Else, see if operator has roots
c too close to the unit circle and should be shrunk.
c-----------------------------------------------------------------------
IF(first)THEN
IF(dpeq(zerom(i),1D0).AND.(.not.onunit))onunit=T
ELSE
IF((zerom(i).le.1.06D0).AND.(.not.shrnkp))shrnkp=T
END IF
END DO
allfix=T
DO ilag=beglag,endlag
IF((.not.Arimaf(ilag)).and.allfix)allfix=F
END DO
c ------------------------------------------------------------------
maon=F
IF(onunit.and.(.not.allfix))THEN
CALL getstr(Oprttl,Oprptr,Noprtl,iopr,tmpttl,ntmpcr)
IF(Lfatal)RETURN
IF(.not.Laumts)WRITE(STDERR,1020)tmpttl(1:ntmpcr)
WRITE(Mt2,1020)tmpttl(1:ntmpcr)
IF(.not.Laumts)WRITE(Mt1,1020)tmpttl(1:ntmpcr)
1020 FORMAT(/,' ERROR: ',a,' polynomial with initial parameters',
& ' is noninvertible',/,' with root(s) on the ',
& 'unit circle. RESPECIFY model with',/,
& ' different initial parameters.',/)
maon=T
maonua=T
END IF
c-----------------------------------------------------------------------
c CODE ADDED BY Brian Monsell Jan. 1998
c Multiply model parameter estimates for this operator by a
c constant, and update the entries in the estimated parameter
c vector.
c ------------------------------------------------------------------
IF(shrnkp.AND.(.not.allfix))THEN
c CALL getstr(Oprttl,Oprptr,Noprtl,iopr,tmpttl,ntmpcr)
c IF(Lfatal)RETURN
c WRITE(STDERR,1021)tmpttl(1:ntmpcr)
c WRITE(Mt2,1021)tmpttl(1:ntmpcr)
c WRITE(Mt1,1021)tmpttl(1:ntmpcr)
c1021 FORMAT(/,' WARNING: ',a,' polynomial from a previous estimation',
c & ' has root(s)',/,' on or near the unit circle.')
DO ilag=beglag,endlag
IF(.not.Arimaf(ilag))THEN
Arimap(ilag)=Arimap(ilag)*(PT9**Arimal(ilag))
Estprm(lagind(ilag))=Arimap(ilag)
END IF
END DO
END IF
c ------------------------------------------------------------------
IF((mains.or.maon).and.(.not.Laumts))THEN
WRITE(Mt1,1030)tmpttl(1:ntmpcr),dotln
1030 FORMAT(' ',a,' Roots',/,' Root',t25,'Real',t31,'Imaginary',
& t44,'Modulus',t53,'Frequency',/,a)
c ------------------------------------------------------------------
DO i=1,degree
WRITE(Mt1,1040)i,zeror(i),zeroi(i),zerom(i),zerof(i)
1040 FORMAT(' Root',i3,t18,4F11.4)
END DO
WRITE(Mt1,1000)dotln
END IF
c-----------------------------------------------------------------------
END DO
END IF
c-----------------------------------------------------------------------
c Compute the roots of initial phi(B)=0
c-----------------------------------------------------------------------
begopr=Mdl(AR-1)
beglag=Opr(begopr-1)
endopr=Mdl(AR)-1
c ------------------------------------------------------------------
IF(endopr.gt.0)THEN
endlag=Opr(endopr)-1
arona=F
c ------------------------------------------------------------------
DO iopr=begopr,endopr
beglag=Opr(iopr-1)
endlag=Opr(iopr)-1
c ------------------------------------------------------------------
factor=Oprfac(iopr)
degree=Arimal(endlag)/factor
coef(1)=-1.0D0
c ------------------------------------------------------------------
CALL setdp(0D0,PORDER,zeror)
CALL setdp(0D0,PORDER,zeroi)
CALL setdp(0D0,PORDER,zerom)
c ------------------------------------------------------------------
CALL setdp(0D0,degree,coef(2))
c ------------------------------------------------------------------
DO ilag=beglag,endlag
coef(Arimal(ilag)/factor+1)=Arimap(ilag)
END DO
c ------------------------------------------------------------------
CALL roots(coef,degree,allinv,zeror,zeroi,zerom,zerof)
c-----------------------------------------------------------------------
c Check stationarity the roots are g(i)=(zeror(i), zeroi(i)),
c i=1,2,...,degree and the complex roots are g(i) and g(i+1).
c If all zeros are stationary do nothing; otherwise print a warning
c message. The program may bomb later if the exact AR is used.
c-----------------------------------------------------------------------
onunit=F
i=0
DO WHILE (i.lt.degree)
i=i+1
IF(dpeq(zerom(i),1D0).and..not.onunit)onunit=T
END DO
IF((.not.allinv).or.onunit)THEN
CALL getstr(Oprttl,Oprptr,Noprtl,iopr,tmpttl,ntmpcr)
IF(Lfatal)RETURN
IF(Lar)THEN
IF(.not.Laumts)WRITE(STDERR,1050)tmpttl(1:ntmpcr)
WRITE(Mt2,1050)tmpttl(1:ntmpcr)
IF(.not.Laumts)WRITE(Mt1,1050)tmpttl(1:ntmpcr)
1050 FORMAT(/,' ERROR: ',a,' polynomial with initial parameters',
& ' is nonstationary',/,' with root(s) on or',
& ' inside the unit circle. RESPECIFY the',/,
& ' model with different initial parameters.',
& /)
arona=T
c ------------------------------------------------------------------
ELSE
IF(.not.(Laumts.or.Lquiet))WRITE(STDERR,1060)tmpttl(1:ntmpcr)
WRITE(Mt2,1060)tmpttl(1:ntmpcr)
IF(.not.Laumts)WRITE(Mt1,1060)tmpttl(1:ntmpcr)
1060 FORMAT(/,' WARNING: ',a,' polynomial with initial parameters',
& ' is nonstationary',
& /,' with root(s) on or inside the unit ',
& 'circle. RESPECIFY the model',
& /,' with different initial parameters.',/)
END IF
c ------------------------------------------------------------------
IF(.not.Laumts)THEN
WRITE(Mt1,1030)tmpttl(1:ntmpcr),dotln
c ------------------------------------------------------------------
DO i=1,degree
WRITE(Mt1,1040)i,zeror(i),zeroi(i),zerom(i),zerof(i)
END DO
WRITE(Mt1,1000)dotln
END IF
c ------------------------------------------------------------------
c IF(.not.arona)arona=T
c ------------------------------------------------------------------
END IF
END DO
END IF
c-----------------------------------------------------------------------
c CODE ADDED BY Brian Monsell Jan. 1998
c ------------------------------------------------------------------
first=.false.
c ------------------------------------------------------------------
IF(mainsa.or.maonua.or.arona)THEN
IF(Laumts)THEN
Laumts=F
ELSE
CALL abend
END IF
END IF
c ------------------------------------------------------------------
1000 FORMAT(a)
c ------------------------------------------------------------------
RETURN
END
|
