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
|
SUBROUTINE tstodf(Trnsrs,Frstry,Nefobs,A,Na,Lsumm,lpr,ldr,lqr,
& lps,lds,lqs,Kstep,Lidotl,Lnoprt,FctOK,redoMD,
& argok)
IMPLICIT NONE
c-----------------------------------------------------------------------
INCLUDE 'stdio.i'
INCLUDE 'notset.prm'
INCLUDE 'srslen.prm'
INCLUDE 'model.prm'
INCLUDE 'model.cmn'
INCLUDE 'arima.cmn'
INCLUDE 'mdldat.cmn'
INCLUDE 'mdldg.cmn'
INCLUDE 'prior.prm'
INCLUDE 'prior.cmn'
INCLUDE 'tbllog.prm'
INCLUDE 'tbllog.cmn'
INCLUDE 'mdltbl.i'
INCLUDE 'svllog.prm'
INCLUDE 'svllog.cmn'
INCLUDE 'mdlsvl.i'
INCLUDE 'units.cmn'
INCLUDE 'error.cmn'
INCLUDE 'extend.cmn'
c-----------------------------------------------------------------------
DOUBLE PRECISION MALIM,ONE,TWO,ZERO,PT5
LOGICAL T,F
PARAMETER(T=.true.,F=.false.,MALIM=0.001D0,ONE=1D0,TWO=2D0,
& ZERO=0D0,PT5=0.05D0)
c-----------------------------------------------------------------------
CHARACTER effttl*(PCOLCR),cmonth*3,ordend*2
DOUBLE PRECISION Trnsrs,sumMA,A,xpxinv,tmp
INTEGER disp,lpr,ldr,lps,lds,lqr,lqs,Frstry,Nefobs,Na,Lsumm,spm1,
& cnote,i,ipos,icol,icol1,igrp,begcol,endcol,regidx,
& nb2,nfix,nchr,nelt,j,Kstep
LOGICAL Lidotl,Lnoprt,FctOK,redoMd,reReDoMd,argok
DIMENSION Trnsrs(PLEN),A(*),cmonth(12),ordend(10),xpxinv(PXPX),
& regidx(PB),tmp(2)
c ------------------------------------------------------------------
DATA cmonth/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep',
& 'Oct','Nov','Dec'/
DATA ordend/'th','st','nd','rd','th','th','th','th','th','th'/
c ------------------------------------------------------------------
DOUBLE PRECISION dpmpar
INTEGER strinx
EXTERNAL dpmpar,strinx
c-----------------------------------------------------------------------
c If model has nonseasonal diffencing and MA, check for
c overdifferencing by seeing if sum of MA parameters is close to 1.
c-----------------------------------------------------------------------
redoMd=F
reReDoMd=F
IF (ldr.gt.0 .and. lqr.gt.0) THEN
IF(Prttab(LAUFNT))write(Mt1,1011)
& ' Checking for nonseasonal overdifferencing.'
disp=lpr+ldr+lps+lds
sumMa=ZERO
do i = disp+1,disp+lqr
sumMa = sumMa + Arimap(i)
end do
c-----------------------------------------------------------------------
c Reduce by one the number of nonseasonal differences, nonseasonal
c MA terms in the model, and add a constant term.
c-----------------------------------------------------------------------
IF (ABS(sumMA-ONE).lt.MALIM) THEN
redoMd=T
IF(Prttab(LAUFNT))write(Mt1,1010)
& ' Reduce order of nonseasonal MA, nonseasonal differencing.'
ldr=ldr-1
lqr=lqr-1
IF(Prttab(LAUFNT))WRITE(Mt1,1020) lpr,ldr,lqr,lps,lds,lqs
icol=strinx(T,Grpttl,Grpptr,1,Ngrptl,'Constant')
IF(icol.eq.0)THEN
IF(Lchkmu)THEN
IF(Prttab(LAUFNT))write(Mt1,1010)' Add constant term.'
CALL adrgef(DNOTST,'Constant','Constant',PRGTCN,F,F)
IF(Lfatal)RETURN
ELSE
cnote=Mt1
IF(Lnoprt.and.(.not.Prttab(LAUFNT)))cnote=0
CALL writln('NOTE: Due to a reduction in the order of '//
& 'regular differencing, a constant term',
& cnote,Mt2,T)
CALL writln(' should be added to the regARIMA model.',
& cnote,Mt2,F)
CALL writln(' Either rerun the spec file with '//
& 'checkmu = yes in the automdl spec, or',
& cnote,Mt2,T)
CALL writln(' add a constant regressor to the '//
& 'regARIMA model via the regression spec.',
& cnote,Mt2,F)
END IF
END IF
END IF
END IF
c-----------------------------------------------------------------------
c If model has seasonal diffencing and MA, check for
c overdifferencing by seeing if sum of seasonal MA parameters is
c close to 1.
c-----------------------------------------------------------------------
IF (lds.gt.0 .and. lqs.gt.0 .and. Lsovdf) THEN
IF(Prttab(LAUFNT))write(Mt1,1011)
& ' Checking for seasonal overdifferencing.'
disp=lpr+ldr+lps+lds+lqr
sumMa=ZERO
do i = disp+1,disp+lqs
sumMa = sumMa + Arimap(i)
end do
c-----------------------------------------------------------------------
c Reduce by one the number of nonseasonal differences, nonseasonal
c MA terms in the model, and add a constant term.
c-----------------------------------------------------------------------
IF (ABS(sumMA-ONE).lt.MALIM) THEN
IF(.not.redoMd)redoMd=T
IF(Prttab(LAUFNT))write(Mt1,1010)
& ' Reduce order of seasonal MA, seasonal differencing.'
lds=lds-1
lqs=lqs-1
IF(Prttab(LAUFNT))WRITE(Mt1,1020) lpr,ldr,lqr,lps,lds,lqs
IF(Prttab(LAUFNT))write(Mt1,1010)' Add seasonal regressors.'
spm1=Sp-1
c ------------------------------------------------------------------
IF(Sp.eq.12)THEN
DO i=1,spm1
effttl=cmonth(i)
nchr=3
CALL adrgef(DNOTST,effttl(1:nchr),'Seasonal',PRGTSE,F,T)
IF(Lfatal)RETURN
END DO
c ------------------------------------------------------------------
ELSE
DO i=1,spm1
ipos=1
CALL itoc(i,effttl,ipos)
IF(Lfatal)RETURN
IF(mod(i,100).ge.11.and.mod(i,100).le.13)THEN
effttl(ipos:ipos+1)='th'
ELSE
effttl(ipos:ipos+1)=ordend(mod(i,10))
END IF
nchr=ipos+1
CALL adrgef(DNOTST,effttl(1:nchr),'Seasonal',PRGTSE,F,T)
IF(Lfatal)RETURN
END DO
END IF
IF(Lfatal)RETURN
END IF
END IF
c-----------------------------------------------------------------------
c re-estimate model.
c-----------------------------------------------------------------------
IF(redoMD)THEN
CALL mdlint()
CALL mdlset(lpr,ldr,lqr,lps,lds,lqs,argok)
IF(.not.Lfatal)
& CALL regvar(Trnsrs,Nobspf,Fctdrp,Nfcst,0,Userx,Bgusrx,Nrusrx,
& Priadj,Reglom,Nrxy,Begxy,Frstry,T,Elong)
IF(Lfatal)RETURN
c-----------------------------------------------------------------------
c If automatic outliers are identified for the model, eliminate the
c outliers from the model (BCM April 2007)
c-----------------------------------------------------------------------
IF(Natotl.gt.0)THEN
CALL clrotl(Nrxy)
IF(.not.Lfatal)
& CALL regvar(Trnsrs,Nobspf,Fctdrp,Nfcst,0,Userx,Bgusrx,
& Nrusrx,Priadj,Reglom,Nrxy,Begxy,Frstry,T,Elong)
IF(Lfatal)RETURN
END IF
c-----------------------------------------------------------------------
CALL rgarma(Lestim,Mxiter,Mxnlit,F,a,na,nefobs,argok)
IF(.not.Lfatal)THEN
CALL prterr(nefobs,T)
IF(.not.Convrg)THEN
WRITE(STDERR,1090)
IF(Prttab(LAUFNT))WRITE(Mt1,1090)
WRITE(Mt2,1090)
CALL abend()
ELSE IF(.not.argok)THEN
CALL abend()
END IF
END IF
IF(Lfatal)RETURN
c-----------------------------------------------------------------------
c Redo automatic outlier identification (BCM April 2007)
c-----------------------------------------------------------------------
IF(Lidotl.and.(.not.Lotmod))THEN
CALL amidot(A,Trnsrs,Frstry,Nefobs,Priadj,Convrg,Fctok,Argok)
IF(Lfatal)RETURN
END IF
c-----------------------------------------------------------------------
c Check if mean term added is significant.
c If not, remove and re-estimate model.
c Added by B C Monsell Aug 2018
c-----------------------------------------------------------------------
IF(Lchkmu)THEN
CALL chkmu(Trnsrs,A,Nefobs,Na,Frstry,kstep,Prttab(LAUFNT))
IF(Lfatal)RETURN
icol1=strinx(T,Grpttl,Grpptr,1,Ngrptl,'Constant')
IF(icol1.eq.0)reReDoMd=T
END IF
c-----------------------------------------------------------------------
c Check if seasonal regressors added are significant.
c If not, remove and re-estimate model.
c Added by B C Monsell Aug 2018
c-----------------------------------------------------------------------
IF(Lsovdf)THEN
icol=strinx(T,Grpttl,Grpptr,1,Ngrptl,'Seasonal')
IF(icol.gt.0)THEN
c ------------------------------------------------------------------
c Generate number of unfixed regressors
c ------------------------------------------------------------------
nb2=Nb
IF(Iregfx.ge.2)THEN
DO j=1,Nb
IF(Regfx(j))nb2=nb2-1
END DO
END IF
c-----------------------------------------------------------------------
c Get the X'X inverse.
c-----------------------------------------------------------------------
IF(nb2.gt.0)THEN
nelt=(nb2+1)*(nb2+2)/2
IF(Var.gt.TWO*dpmpar(1))THEN
CALL copy(Chlxpx,nelt,1,xpxinv)
CALL dppdi(xpxinv,nb2,tmp,1)
END IF
END IF
c-----------------------------------------------------------------------
c set up regidx variable
c-----------------------------------------------------------------------
nfix=0
DO igrp=1,Ngrp
begcol=Grp(igrp-1)
endcol=Grp(igrp)-1
DO icol=begcol,endcol
IF(Regfx(icol))THEN
nfix=nfix+1
regidx(icol)=NOTSET
ELSE
regidx(icol)=icol-nfix
END IF
END DO
END DO
c-----------------------------------------------------------------------
c Generate F-test for seasonal F-test
c-----------------------------------------------------------------------
CALL sftest(Xpxinv,Regidx,Prttab(LAUFNT),F,F,F)
c-----------------------------------------------------------------------
c generate p-value limit,
c remove seasonal regressors if not significant
c-----------------------------------------------------------------------
IF(Sfpv.gt.PT5)THEN
igrp=strinx(T,Grpttl,Grpptr,1,Ngrptl,'Seasonal')
begcol=Grp(igrp-1)
endcol=Grp(igrp)-1
CALL dlrgef(begcol,Nrxy,endcol-begcol+1)
IF(Prttab(LAUFNT))
& write(Mt1,1010)'Seasonal regressors removed from model'
END IF
icol1=strinx(T,Grpttl,Grpptr,1,Ngrptl,'Seasonal')
IF(icol1.eq.0)reReDoMd=T
END IF
END IF
c-----------------------------------------------------------------------
c Re-estimate model if
c-----------------------------------------------------------------------
IF(reReDoMd)THEN
CALL rgarma(Lestim,Mxiter,Mxnlit,F,a,na,nefobs,argok)
IF(.not.Lfatal)THEN
CALL prterr(nefobs,T)
IF(.not.Convrg)THEN
WRITE(STDERR,1090)
IF(Prttab(LAUFNT))WRITE(Mt1,1090)
WRITE(Mt2,1090)
CALL abend()
ELSE IF(.not.argok)THEN
CALL abend()
END IF
END IF
IF(Lfatal)RETURN
END IF
c-----------------------------------------------------------------------
ELSE
IF(Prttab(LAUFNT))WRITE(Mt1,1040)MALIM
END IF
c-----------------------------------------------------------------------
1010 FORMAT(' ',a)
1011 FORMAT(/,' ',a)
1020 FORMAT(' ',2(' (',i2,',',i2,',',i2,')'))
1090 FORMAT(/,' ERROR: Estimation failed to converge during the ',
& 'automatic model',
& /,' identification procedure.')
1040 FORMAT(' Nonseasonal MA not within ',f6.3,' of 1.0 - model ',
& 'passes test.')
c-----------------------------------------------------------------------
RETURN
END
|
