1c
2c     from netlib/a/stl: no authorship nor copyright claim in the source;
3c     presumably by the authors of
4c
5c     R.B. Cleveland, W.S.Cleveland, J.E. McRae, and I. Terpenning,
6c     STL: A Seasonal-Trend Decomposition Procedure Based on Loess,
7c     Statistics Research Report, AT&T Bell Laboratories.
8c
9c     Converted to double precision by B.D. Ripley 1999.
10c     Indented, goto labels renamed, many goto's replaced by `if then {else}'
11c     (using Emacs), many more comments;  by M.Maechler 2001-02.
12c
13      subroutine stl(y,n,np,ns,nt,nl, isdeg,itdeg,ildeg,
14     &     nsjump,ntjump,nljump, ni,no, rw,season,trend,work)
15
16c     implicit none
17c Arg
18      integer n, np, ns,nt,nl, isdeg,itdeg,ildeg, nsjump,ntjump,nljump,
19     &     ni, no
20c       n                   : length(y)
21c       ns, nt, nl          : spans        for `s', `t' and `l' smoother
22c       isdeg, itdeg, ildeg : local degree for `s', `t' and `l' smoother
23c       nsjump,ntjump,nljump: ........     for `s', `t' and `l' smoother
24c       ni, no              : number of inner and outer (robust) iterations
25
26      double precision y(n), rw(n), season(n), trend(n),
27     &     work(n+2*np,5)
28c Var
29      integer i,k, newns, newnt, newnl, newnp
30      logical userw
31
32      userw = .false.
33      do 1 i = 1,n
34         trend(i) = 0.d0
35 1    continue
36c the three spans must be at least three and odd:
37      newns = max0(3,ns)
38      newnt = max0(3,nt)
39      newnl = max0(3,nl)
40      if(mod(newns,2) .eq. 0) newns = newns + 1
41      if(mod(newnt,2) .eq. 0) newnt = newnt + 1
42      if(mod(newnl,2) .eq. 0) newnl = newnl + 1
43c periodicity at least 2:
44      newnp = max0(2,np)
45
46      k = 0
47c --- outer loop -- robustnes iterations
48 100  continue
49      call stlstp(y,n, newnp,newns,newnt,newnl, isdeg,itdeg,ildeg,
50     &     nsjump,ntjump,nljump, ni,userw,rw,season, trend, work)
51      k = k+1
52      if(k .gt. no) goto 10
53
54      do 3 i = 1,n
55         work(i,1) = trend(i)+season(i)
56 3    continue
57      call stlrwt(y,n,work(1,1),rw)
58      userw = .true.
59      goto 100
60c --- end Loop
61 10   continue
62
63c     robustness weights when there were no robustness iterations:
64      if(no .le. 0) then
65         do 15 i = 1,n
66            rw(i) = 1.d0
67 15      continue
68      endif
69      return
70      end
71
72      subroutine stless(y,n,len,ideg,njump,userw,rw,ys,res)
73
74c     implicit none
75c Arg
76      integer n, len, ideg, njump
77      double precision y(n), rw(n), ys(n), res(n)
78c Var
79      integer newnj, nleft, nright, nsh, k, i, j
80      double precision delta
81      logical ok, userw
82
83      if(n .lt. 2) then
84         ys(1) = y(1)
85         return
86      endif
87
88      newnj = min0(njump, n-1)
89      if(len .ge. n) then
90         nleft = 1
91         nright = n
92         do 20 i = 1,n,newnj
93            call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,
94     &           userw,rw,ok)
95            if(.not. ok) ys(i) = y(i)
96 20      continue
97
98      else
99
100         if(newnj .eq. 1) then
101            nsh = (len+1)/2
102            nleft = 1
103            nright = len
104            do 30 i = 1,n
105               if(i .gt. nsh  .and.  nright .ne. n) then
106                  nleft = nleft+1
107                  nright = nright+1
108               endif
109               call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,
110     &              userw,rw,ok)
111               if(.not. ok) ys(i) = y(i)
112 30         continue
113         else
114            nsh = (len+1)/2
115            do 40 i = 1,n,newnj
116               if(i .lt. nsh) then
117                  nleft = 1
118                  nright = len
119               else if(i .ge. n-nsh+1) then
120                  nleft = n-len+1
121                  nright = n
122               else
123                  nleft = i-nsh+1
124                  nright = len+i-nsh
125               endif
126
127               call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,
128     &              userw,rw,ok)
129               if(.not. ok) ys(i) = y(i)
130 40         continue
131
132         endif
133
134      endif
135
136      if(newnj .ne. 1) then
137         do 45 i = 1,n-newnj,newnj
138            delta = (ys(i+newnj)-ys(i))/dble(newnj)
139            do 47 j = i+1,i+newnj-1
140               ys(j) = ys(i)+delta*dble(j-i)
141 47         continue
142 45      continue
143         k = ((n-1)/newnj)*newnj+1
144
145         if(k .ne. n) then
146            call stlest(y,n,len,ideg,dble(n),ys(n),nleft,nright,res,
147     &           userw,rw,ok)
148            if(.not. ok) ys(n) = y(n)
149
150            if(k .ne. n-1) then
151               delta = (ys(n)-ys(k))/dble(n-k)
152               do 55 j = k+1,n-1
153                  ys(j) = ys(k)+delta*dble(j-k)
154 55            continue
155            endif
156         endif
157      endif
158      return
159      end
160
161      subroutine stlest(y,n,len,ideg,xs,ys,nleft,nright,w,
162     &     userw,rw,ok)
163
164c     implicit none
165c Arg
166      integer n, len, ideg, nleft, nright
167      double precision y(n), w(n), rw(n), xs, ys
168      logical userw,ok
169c Var
170      double precision range, h, h1, h9, a, b, c, r
171      integer j
172
173      range = dble(n)-dble(1)
174      h = max(xs - dble(nleft), dble(nright) - xs)
175      if(len .gt. n) h = h + dble((len-n)/2)
176      h9 = 0.999d0*h
177      h1 = 0.001d0*h
178      a = 0.d0
179      do 60 j = nleft,nright
180         r = abs(dble(j)-xs)
181         if(r .le. h9) then
182            if(r .le. h1) then
183               w(j) = 1.d0
184            else
185               w(j) = (1.d0 - (r/h)**3)**3
186            endif
187            if(userw) w(j) = rw(j)*w(j)
188            a = a+w(j)
189         else
190            w(j) = 0.d0
191         endif
192 60   continue
193
194      if(a .le. 0.d0) then
195         ok = .false.
196      else
197         ok = .true.
198         do 69 j = nleft,nright
199            w(j) = w(j)/a
200 69      continue
201         if((h .gt. 0.d0) .and. (ideg .gt. 0)) then
202            a = 0.d0
203            do 73 j = nleft,nright
204               a = a+w(j)*dble(j)
205 73         continue
206            b = xs-a
207            c = 0.d0
208            do 75 j = nleft,nright
209               c = c+w(j)*(dble(j)-a)**2
210 75         continue
211            if(sqrt(c) .gt. 0.001d0*range) then
212               b = b/c
213               do 79 j = nleft,nright
214                  w(j) = w(j)*(b*(dble(j)-a)+1.0d0)
215 79            continue
216            endif
217         endif
218         ys = 0.d0
219         do 81 j = nleft,nright
220            ys = ys+w(j)*y(j)
221 81      continue
222      endif
223
224      return
225      end
226
227      subroutine stlfts(x,n,np,trend,work)
228      integer n, np
229      double precision x(n), trend(n), work(n)
230
231      call stlma(x,    n,      np, trend)
232      call stlma(trend,n-np+1, np, work)
233      call stlma(work, n-2*np+2,3, trend)
234      return
235      end
236
237
238      subroutine stlma(x, n, len, ave)
239
240c Moving Average (aka "running mean")
241c ave(i) := mean(x{j}, j = max(1,i-k),..., min(n, i+k))
242c           for i = 1,2,..,n
243
244c     implicit none
245c Arg
246      integer n, len
247      double precision x(n), ave(n)
248c Var
249      double precision flen, v
250      integer i, j, k, m, newn
251      newn = n-len+1
252      flen = dble(len)
253      v = 0.d0
254      do 3 i = 1,len
255         v = v+x(i)
256 3    continue
257      ave(1) = v/flen
258      if(newn .gt. 1) then
259         k = len
260         m = 0
261         do 7 j = 2, newn
262            k = k+1
263            m = m+1
264            v = v-x(m)+x(k)
265            ave(j) = v/flen
266 7       continue
267      endif
268      return
269      end
270
271
272      subroutine stlstp(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,
273     &     ntjump,nljump,ni,userw,rw,season,trend,work)
274
275c     implicit none
276c Arg
277      integer n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,ni
278      logical userw
279      double precision y(n),rw(n),season(n),trend(n),work(n+2*np,5)
280c Var
281      integer i,j
282
283      do 80 j = 1,ni
284         do 1 i = 1,n
285            work(i,1) = y(i)-trend(i)
286 1       continue
287         call stlss(work(1,1),n,np,ns,isdeg,nsjump,userw,rw,work(1,2),
288     &        work(1,3),work(1,4),work(1,5),season)
289         call stlfts(work(1,2),n+2*np,np,work(1,3),work(1,1))
290         call stless(work(1,3),n,nl,ildeg,nljump,.false.,work(1,4),
291     &        work(1,1),work(1,5))
292         do 3 i = 1,n
293            season(i) = work(np+i,2)-work(i,1)
294 3       continue
295         do 5 i = 1,n
296            work(i,1) = y(i)-season(i)
297 5       continue
298         call stless(work(1,1),n,nt,itdeg,ntjump,userw,rw,trend,
299     &        work(1,3))
300 80   continue
301      return
302      end
303
304      subroutine stlrwt(y,n,fit,rw)
305c Robustness Weights
306c       rw_i := B( |y_i - fit_i| / (6 M) ),   i = 1,2,...,n
307c               where B(u) = (1 - u^2)^2  * 1[|u| < 1]   {Tukey's biweight}
308c               and   M := median{ |y_i - fit_i| }
309c     implicit none
310c Arg
311      integer n
312      double precision y(n), fit(n), rw(n)
313c Var
314      integer mid(2), i
315      double precision cmad, c9, c1, r
316
317      do 7 i = 1,n
318         rw(i) = abs(y(i)-fit(i))
319 7    continue
320      mid(1) = n/2+1
321      mid(2) = n-mid(1)+1
322      call psort(rw,n,mid,2)
323      cmad = 3.0d0*(rw(mid(1))+rw(mid(2)))
324c     = 6 * MAD
325      c9 = 0.999d0*cmad
326      c1 = 0.001d0*cmad
327      do 10 i = 1,n
328         r = abs(y(i)-fit(i))
329         if(r .le. c1) then
330            rw(i) = 1.d0
331         else if(r .le. c9) then
332            rw(i) = (1.d0 - (r/cmad)**2)**2
333         else
334            rw(i) = 0.d0
335         endif
336 10   continue
337      return
338      end
339
340      subroutine stlss(y,n,np,ns,isdeg,nsjump,userw,rw,season,
341     &     work1,work2,work3,work4)
342c
343c       called by stlstp() at the beginning of each (inner) iteration
344c
345c     implicit none
346c Arg
347      integer n, np, ns, isdeg, nsjump
348      double precision y(n), rw(n), season(n+2*np),
349     &     work1(n), work2(n), work3(n), work4(n)
350      logical userw
351c Var
352      integer nright, nleft, i, j, k, m
353      logical ok
354      double precision xs
355
356      if(np .lt. 1) return
357
358      do 200 j = 1, np
359         k = (n-j)/np+1
360         do 10 i = 1,k
361            work1(i) = y((i-1)*np+j)
362 10      continue
363         if(userw) then
364            do 12 i = 1,k
365               work3(i) = rw((i-1)*np+j)
366 12         continue
367         endif
368         call stless(work1,k,ns,isdeg,nsjump,userw,work3,work2(2),work4)
369         xs = 0
370         nright = min0(ns,k)
371         call stlest(work1,k,ns,isdeg,xs,work2(1),1,nright,work4,
372     &        userw,work3,ok)
373         if(.not. ok) work2(1) = work2(2)
374         xs = k+1
375         nleft = max0(1,k-ns+1)
376         call stlest(work1,k,ns,isdeg,xs,work2(k+2),nleft,k,work4,
377     &        userw,work3,ok)
378         if(.not. ok) work2(k+2) = work2(k+1)
379         do 18 m = 1,k+2
380            season((m-1)*np+j) = work2(m)
381 18      continue
382
383 200  continue
384
385      return
386      end
387
388
389c  STL E_Z_ : "Easy" user interface  -- not called from R
390
391      subroutine stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw,
392     &     season, trend, work)
393
394c     implicit none
395c Arg
396      integer n, np, ns, isdeg, itdeg, no
397      logical robust
398      double precision y(n), rw(n), season(n), trend(n), work(n+2*np,7)
399c Var
400      integer i, j, ildeg, nt, nl, ni, nsjump, ntjump, nljump,
401     &     newns, newnp
402      double precision maxs, mins, maxt, mint, maxds, maxdt, difs, dift
403
404      ildeg = itdeg
405      newns = max0(3,ns)
406      if(mod(newns,2) .eq. 0) newns = newns+1
407      newnp = max0(2,np)
408      nt = int((1.5d0*newnp)/(1.d0 - 1.5d0/newns) + 0.5d0)
409      nt = max0(3,nt)
410      if(mod(nt,2) .eq. 0) nt = nt+1
411      nl = newnp
412      if(mod(nl,2) .eq. 0) nl = nl+1
413
414      if(robust) then
415         ni = 1
416      else
417         ni = 2
418      endif
419
420      nsjump = max0(1,int(float(newns)/10 + 0.9))
421      ntjump = max0(1,int(float(nt)/10 + 0.9))
422      nljump = max0(1,int(float(nl)/10 + 0.9))
423      do 2 i = 1,n
424         trend(i) = 0.d0
425 2    continue
426      call stlstp(y,n,newnp,newns,nt,nl,isdeg,itdeg,ildeg,nsjump,
427     &     ntjump,nljump,ni,.false.,rw,season,trend,work)
428
429      no = 0
430      if(robust) then
431         j=1
432C        Loop  --- 15 robustness iterations
433 100     if(j .le. 15) then
434            do 35 i = 1,n
435               work(i,6) = season(i)
436               work(i,7) = trend(i)
437               work(i,1) = trend(i)+season(i)
438 35        continue
439            call stlrwt(y,n,work(1,1),rw)
440            call stlstp(y, n, newnp, newns, nt,nl, isdeg,itdeg,ildeg,
441     &           nsjump,ntjump,nljump, ni, .true.,
442     &           rw, season, trend, work)
443            no = no+1
444            maxs = work(1,6)
445            mins = work(1,6)
446            maxt = work(1,7)
447            mint = work(1,7)
448            maxds = abs(work(1,6) - season(1))
449            maxdt = abs(work(1,7) - trend(1))
450            do   137 i = 2,n
451               if(maxs .lt. work(i,6)) maxs = work(i,6)
452               if(maxt .lt. work(i,7)) maxt = work(i,7)
453               if(mins .gt. work(i,6)) mins = work(i,6)
454               if(mint .gt. work(i,7)) mint = work(i,7)
455               difs = abs(work(i,6) - season(i))
456               dift = abs(work(i,7) - trend(i))
457               if(maxds .lt. difs) maxds = difs
458               if(maxdt .lt. dift) maxdt = dift
459 137        continue
460            if((maxds/(maxs-mins) .lt. 0.01d0) .and.
461     &         (maxdt/(maxt-mint) .lt. 0.01d0))   goto 300
462            continue
463            j=j+1
464            goto 100
465         endif
466C        end Loop
467 300     continue
468
469      else
470c       .not. robust
471
472         do 150 i = 1,n
473            rw(i) = 1.0d0
474 150     continue
475      endif
476
477      return
478      end
479
480      subroutine psort(a,n,ind,ni)
481c
482c Partial Sorting ; used for Median (MAD) computation only
483c
484c     implicit none
485c Arg
486      integer n,ni
487      double precision a(n)
488      integer ind(ni)
489c Var
490      integer indu(16),indl(16),iu(16),il(16),p,jl,ju,i,j,m,k,ij,l
491      double precision t,tt
492
493      if(n .lt. 0 .or. ni .lt. 0) return
494
495      if(n .lt. 2 .or. ni .eq. 0) return
496
497      jl = 1
498      ju = ni
499      indl(1) = 1
500      indu(1) = ni
501      i = 1
502      j = n
503      m = 1
504
505c Outer Loop
506 161  continue
507      if(i .lt. j) go to 10
508
509c  _Loop_
510 166  continue
511      m = m-1
512      if(m .eq. 0) return
513      i = il(m)
514      j = iu(m)
515      jl = indl(m)
516      ju = indu(m)
517      if(.not.(jl .le. ju))  goto 166
518
519c     while (j - i > 10)
520 173  if(.not.(j-i .gt. 10)) goto 174
521
522 10   k = i
523      ij = (i+j)/2
524      t = a(ij)
525      if(a(i) .gt. t) then
526         a(ij) = a(i)
527         a(i) = t
528         t = a(ij)
529      endif
530      l = j
531      if(a(j) .lt. t) then
532         a(ij) = a(j)
533         a(j) = t
534         t = a(ij)
535         if(a(i) .gt. t) then
536            a(ij) = a(i)
537            a(i) = t
538            t = a(ij)
539         endif
540      endif
541
542 181  continue
543      l = l-1
544      if(a(l) .le. t)then
545         tt = a(l)
546 186     continue
547         k = k+1
548         if(.not.(a(k) .ge. t)) goto 186
549
550         if(k .gt. l) goto 183
551
552         a(l) = a(k)
553         a(k) = tt
554      endif
555      goto   181
556
557 183  continue
558      indl(m) = jl
559      indu(m) = ju
560      p = m
561      m = m+1
562      if(l-i .le. j-k) then
563         il(p) = k
564         iu(p) = j
565         j = l
566
567 193     continue
568         if(jl .gt. ju) goto 166
569         if(ind(ju) .gt. j) then
570            ju = ju-1
571            goto 193
572         endif
573         indl(p) = ju+1
574      else
575         il(p) = i
576         iu(p) = l
577         i = k
578
579 200     continue
580         if(jl .gt. ju) goto 166
581         if(ind(jl) .lt. i) then
582            jl = jl+1
583            goto 200
584         endif
585         indu(p) = jl-1
586      endif
587
588      goto   173
589c     end while
590 174  continue
591
592      if(i .ne. 1) then
593         i = i-1
594 209     continue
595         i = i+1
596         if(i .eq. j) goto 166
597         t = a(i+1)
598         if(a(i) .gt. t) then
599            k = i
600c           repeat
601 216        continue
602            a(k+1) = a(k)
603            k = k-1
604            if(.not.(t .ge. a(k))) goto 216
605c           until  t >= a(k)
606            a(k+1) = t
607         endif
608         goto 209
609
610      endif
611
612      goto 161
613c End Outer Loop
614
615      end
616