| c |
| c from netlib/a/stl: no authorship nor copyright claim in the source; |
| c presumably by the authors of |
| c |
| c R.B. Cleveland, W.S.Cleveland, J.E. McRae, and I. Terpenning, |
| c STL: A Seasonal-Trend Decomposition Procedure Based on Loess, |
| c Statistics Research Report, AT&T Bell Laboratories. |
| c |
| c Converted to double precision by B.D. Ripley 1999. |
| c Indented, goto labels renamed, many goto's replaced by `if then {else}' |
| c (using Emacs), many more comments; by M.Maechler 2001-02. |
| c |
| subroutine stl(y,n,np,ns,nt,nl, isdeg,itdeg,ildeg, |
| & nsjump,ntjump,nljump, ni,no, rw,season,trend,work) |
| |
| c implicit none |
| c Arg |
| integer n, np, ns,nt,nl, isdeg,itdeg,ildeg, nsjump,ntjump,nljump, |
| & ni, no |
| c n : length(y) |
| c ns, nt, nl : spans for `s', `t' and `l' smoother |
| c isdeg, itdeg, ildeg : local degree for `s', `t' and `l' smoother |
| c nsjump,ntjump,nljump: ........ for `s', `t' and `l' smoother |
| c ni, no : number of inner and outer (robust) iterations |
| |
| double precision y(n), rw(n), season(n), trend(n), |
| & work(n+2*np,5) |
| c Var |
| integer i,k, newns, newnt, newnl, newnp |
| logical userw |
| |
| userw = .false. |
| do 1 i = 1,n |
| trend(i) = 0.d0 |
| 1 continue |
| c the three spans must be at least three and odd: |
| newns = max0(3,ns) |
| newnt = max0(3,nt) |
| newnl = max0(3,nl) |
| if(mod(newns,2) .eq. 0) newns = newns + 1 |
| if(mod(newnt,2) .eq. 0) newnt = newnt + 1 |
| if(mod(newnl,2) .eq. 0) newnl = newnl + 1 |
| c periodicity at least 2: |
| newnp = max0(2,np) |
| |
| k = 0 |
| c --- outer loop -- robustnes iterations |
| 100 continue |
| call stlstp(y,n, newnp,newns,newnt,newnl, isdeg,itdeg,ildeg, |
| & nsjump,ntjump,nljump, ni,userw,rw,season, trend, work) |
| k = k+1 |
| if(k .gt. no) goto 10 |
| |
| do 3 i = 1,n |
| work(i,1) = trend(i)+season(i) |
| 3 continue |
| call stlrwt(y,n,work(1,1),rw) |
| userw = .true. |
| goto 100 |
| c --- end Loop |
| 10 continue |
| |
| c robustness weights when there were no robustness iterations: |
| if(no .le. 0) then |
| do 15 i = 1,n |
| rw(i) = 1.d0 |
| 15 continue |
| endif |
| return |
| end |
| |
| subroutine stless(y,n,len,ideg,njump,userw,rw,ys,res) |
| |
| c implicit none |
| c Arg |
| integer n, len, ideg, njump |
| double precision y(n), rw(n), ys(n), res(n) |
| c Var |
| integer newnj, nleft, nright, nsh, k, i, j |
| double precision delta |
| logical ok, userw |
| |
| if(n .lt. 2) then |
| ys(1) = y(1) |
| return |
| endif |
| |
| newnj = min0(njump, n-1) |
| if(len .ge. n) then |
| nleft = 1 |
| nright = n |
| do 20 i = 1,n,newnj |
| call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res, |
| & userw,rw,ok) |
| if(.not. ok) ys(i) = y(i) |
| 20 continue |
| |
| else |
| |
| if(newnj .eq. 1) then |
| nsh = (len+1)/2 |
| nleft = 1 |
| nright = len |
| do 30 i = 1,n |
| if(i .gt. nsh .and. nright .ne. n) then |
| nleft = nleft+1 |
| nright = nright+1 |
| endif |
| call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res, |
| & userw,rw,ok) |
| if(.not. ok) ys(i) = y(i) |
| 30 continue |
| else |
| nsh = (len+1)/2 |
| do 40 i = 1,n,newnj |
| if(i .lt. nsh) then |
| nleft = 1 |
| nright = len |
| else if(i .ge. n-nsh+1) then |
| nleft = n-len+1 |
| nright = n |
| else |
| nleft = i-nsh+1 |
| nright = len+i-nsh |
| endif |
| |
| call stlest(y,n,len,ideg,dble(i),ys(i),nleft,nright,res, |
| & userw,rw,ok) |
| if(.not. ok) ys(i) = y(i) |
| 40 continue |
| |
| endif |
| |
| endif |
| |
| if(newnj .ne. 1) then |
| do 45 i = 1,n-newnj,newnj |
| delta = (ys(i+newnj)-ys(i))/dble(newnj) |
| do 47 j = i+1,i+newnj-1 |
| ys(j) = ys(i)+delta*dble(j-i) |
| 47 continue |
| 45 continue |
| k = ((n-1)/newnj)*newnj+1 |
| |
| if(k .ne. n) then |
| call stlest(y,n,len,ideg,dble(n),ys(n),nleft,nright,res, |
| & userw,rw,ok) |
| if(.not. ok) ys(n) = y(n) |
| |
| if(k .ne. n-1) then |
| delta = (ys(n)-ys(k))/dble(n-k) |
| do 55 j = k+1,n-1 |
| ys(j) = ys(k)+delta*dble(j-k) |
| 55 continue |
| endif |
| endif |
| endif |
| return |
| end |
| |
| subroutine stlest(y,n,len,ideg,xs,ys,nleft,nright,w, |
| & userw,rw,ok) |
| |
| c implicit none |
| c Arg |
| integer n, len, ideg, nleft, nright |
| double precision y(n), w(n), rw(n), xs, ys |
| logical userw,ok |
| c Var |
| double precision range, h, h1, h9, a, b, c, r |
| integer j |
| |
| range = dble(n)-dble(1) |
| h = max(xs - dble(nleft), dble(nright) - xs) |
| if(len .gt. n) h = h + dble((len-n)/2) |
| h9 = 0.999d0*h |
| h1 = 0.001d0*h |
| a = 0.d0 |
| do 60 j = nleft,nright |
| r = abs(dble(j)-xs) |
| if(r .le. h9) then |
| if(r .le. h1) then |
| w(j) = 1.d0 |
| else |
| w(j) = (1.d0 - (r/h)**3)**3 |
| endif |
| if(userw) w(j) = rw(j)*w(j) |
| a = a+w(j) |
| else |
| w(j) = 0.d0 |
| endif |
| 60 continue |
| |
| if(a .le. 0.d0) then |
| ok = .false. |
| else |
| ok = .true. |
| do 69 j = nleft,nright |
| w(j) = w(j)/a |
| 69 continue |
| if((h .gt. 0.d0) .and. (ideg .gt. 0)) then |
| a = 0.d0 |
| do 73 j = nleft,nright |
| a = a+w(j)*dble(j) |
| 73 continue |
| b = xs-a |
| c = 0.d0 |
| do 75 j = nleft,nright |
| c = c+w(j)*(dble(j)-a)**2 |
| 75 continue |
| if(sqrt(c) .gt. 0.001d0*range) then |
| b = b/c |
| do 79 j = nleft,nright |
| w(j) = w(j)*(b*(dble(j)-a)+1.0d0) |
| 79 continue |
| endif |
| endif |
| ys = 0.d0 |
| do 81 j = nleft,nright |
| ys = ys+w(j)*y(j) |
| 81 continue |
| endif |
| |
| return |
| end |
| |
| subroutine stlfts(x,n,np,trend,work) |
| integer n, np |
| double precision x(n), trend(n), work(n) |
| |
| call stlma(x, n, np, trend) |
| call stlma(trend,n-np+1, np, work) |
| call stlma(work, n-2*np+2,3, trend) |
| return |
| end |
| |
| |
| subroutine stlma(x, n, len, ave) |
| |
| c Moving Average (aka "running mean") |
| c ave(i) := mean(x{j}, j = max(1,i-k),..., min(n, i+k)) |
| c for i = 1,2,..,n |
| |
| c implicit none |
| c Arg |
| integer n, len |
| double precision x(n), ave(n) |
| c Var |
| double precision flen, v |
| integer i, j, k, m, newn |
| newn = n-len+1 |
| flen = dble(len) |
| v = 0.d0 |
| do 3 i = 1,len |
| v = v+x(i) |
| 3 continue |
| ave(1) = v/flen |
| if(newn .gt. 1) then |
| k = len |
| m = 0 |
| do 7 j = 2, newn |
| k = k+1 |
| m = m+1 |
| v = v-x(m)+x(k) |
| ave(j) = v/flen |
| 7 continue |
| endif |
| return |
| end |
| |
| |
| subroutine stlstp(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump, |
| & ntjump,nljump,ni,userw,rw,season,trend,work) |
| |
| c implicit none |
| c Arg |
| integer n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,ni |
| logical userw |
| double precision y(n),rw(n),season(n),trend(n),work(n+2*np,5) |
| c Var |
| integer i,j |
| |
| do 80 j = 1,ni |
| do 1 i = 1,n |
| work(i,1) = y(i)-trend(i) |
| 1 continue |
| call stlss(work(1,1),n,np,ns,isdeg,nsjump,userw,rw,work(1,2), |
| & work(1,3),work(1,4),work(1,5),season) |
| call stlfts(work(1,2),n+2*np,np,work(1,3),work(1,1)) |
| call stless(work(1,3),n,nl,ildeg,nljump,.false.,work(1,4), |
| & work(1,1),work(1,5)) |
| do 3 i = 1,n |
| season(i) = work(np+i,2)-work(i,1) |
| 3 continue |
| do 5 i = 1,n |
| work(i,1) = y(i)-season(i) |
| 5 continue |
| call stless(work(1,1),n,nt,itdeg,ntjump,userw,rw,trend, |
| & work(1,3)) |
| 80 continue |
| return |
| end |
| |
| subroutine stlrwt(y,n,fit,rw) |
| c Robustness Weights |
| c rw_i := B( |y_i - fit_i| / (6 M) ), i = 1,2,...,n |
| c where B(u) = (1 - u^2)^2 * 1[|u| < 1] {Tukey's biweight} |
| c and M := median{ |y_i - fit_i| } |
| c implicit none |
| c Arg |
| integer n |
| double precision y(n), fit(n), rw(n) |
| c Var |
| integer mid(2), i |
| double precision cmad, c9, c1, r |
| |
| do 7 i = 1,n |
| rw(i) = abs(y(i)-fit(i)) |
| 7 continue |
| mid(1) = n/2+1 |
| mid(2) = n-mid(1)+1 |
| call psort(rw,n,mid,2) |
| cmad = 3.0d0*(rw(mid(1))+rw(mid(2))) |
| c = 6 * MAD |
| c9 = 0.999d0*cmad |
| c1 = 0.001d0*cmad |
| do 10 i = 1,n |
| r = abs(y(i)-fit(i)) |
| if(r .le. c1) then |
| rw(i) = 1.d0 |
| else if(r .le. c9) then |
| rw(i) = (1.d0 - (r/cmad)**2)**2 |
| else |
| rw(i) = 0.d0 |
| endif |
| 10 continue |
| return |
| end |
| |
| subroutine stlss(y,n,np,ns,isdeg,nsjump,userw,rw,season, |
| & work1,work2,work3,work4) |
| c |
| c called by stlstp() at the beginning of each (inner) iteration |
| c |
| c implicit none |
| c Arg |
| integer n, np, ns, isdeg, nsjump |
| double precision y(n), rw(n), season(n+2*np), |
| & work1(n), work2(n), work3(n), work4(n) |
| logical userw |
| c Var |
| integer nright, nleft, i, j, k, m |
| logical ok |
| double precision xs |
| |
| if(np .lt. 1) return |
| |
| do 200 j = 1, np |
| k = (n-j)/np+1 |
| do 10 i = 1,k |
| work1(i) = y((i-1)*np+j) |
| 10 continue |
| if(userw) then |
| do 12 i = 1,k |
| work3(i) = rw((i-1)*np+j) |
| 12 continue |
| endif |
| call stless(work1,k,ns,isdeg,nsjump,userw,work3,work2(2),work4) |
| xs = 0 |
| nright = min0(ns,k) |
| call stlest(work1,k,ns,isdeg,xs,work2(1),1,nright,work4, |
| & userw,work3,ok) |
| if(.not. ok) work2(1) = work2(2) |
| xs = k+1 |
| nleft = max0(1,k-ns+1) |
| call stlest(work1,k,ns,isdeg,xs,work2(k+2),nleft,k,work4, |
| & userw,work3,ok) |
| if(.not. ok) work2(k+2) = work2(k+1) |
| do 18 m = 1,k+2 |
| season((m-1)*np+j) = work2(m) |
| 18 continue |
| |
| 200 continue |
| |
| return |
| end |
| |
| |
| c STL E_Z_ : "Easy" user interface -- not called from R |
| |
| subroutine stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, |
| & season, trend, work) |
| |
| c implicit none |
| c Arg |
| integer n, np, ns, isdeg, itdeg, no |
| logical robust |
| double precision y(n), rw(n), season(n), trend(n), work(n+2*np,7) |
| c Var |
| integer i, j, ildeg, nt, nl, ni, nsjump, ntjump, nljump, |
| & newns, newnp |
| double precision maxs, mins, maxt, mint, maxds, maxdt, difs, dift |
| |
| ildeg = itdeg |
| newns = max0(3,ns) |
| if(mod(newns,2) .eq. 0) newns = newns+1 |
| newnp = max0(2,np) |
| nt = int((1.5d0*newnp)/(1.d0 - 1.5d0/newns) + 0.5d0) |
| nt = max0(3,nt) |
| if(mod(nt,2) .eq. 0) nt = nt+1 |
| nl = newnp |
| if(mod(nl,2) .eq. 0) nl = nl+1 |
| |
| if(robust) then |
| ni = 1 |
| else |
| ni = 2 |
| endif |
| |
| nsjump = max0(1,int(float(newns)/10 + 0.9)) |
| ntjump = max0(1,int(float(nt)/10 + 0.9)) |
| nljump = max0(1,int(float(nl)/10 + 0.9)) |
| do 2 i = 1,n |
| trend(i) = 0.d0 |
| 2 continue |
| call stlstp(y,n,newnp,newns,nt,nl,isdeg,itdeg,ildeg,nsjump, |
| & ntjump,nljump,ni,.false.,rw,season,trend,work) |
| |
| no = 0 |
| if(robust) then |
| j=1 |
| C Loop --- 15 robustness iterations |
| 100 if(j .le. 15) then |
| do 35 i = 1,n |
| work(i,6) = season(i) |
| work(i,7) = trend(i) |
| work(i,1) = trend(i)+season(i) |
| 35 continue |
| call stlrwt(y,n,work(1,1),rw) |
| call stlstp(y, n, newnp, newns, nt,nl, isdeg,itdeg,ildeg, |
| & nsjump,ntjump,nljump, ni, .true., |
| & rw, season, trend, work) |
| no = no+1 |
| maxs = work(1,6) |
| mins = work(1,6) |
| maxt = work(1,7) |
| mint = work(1,7) |
| maxds = abs(work(1,6) - season(1)) |
| maxdt = abs(work(1,7) - trend(1)) |
| do 137 i = 2,n |
| if(maxs .lt. work(i,6)) maxs = work(i,6) |
| if(maxt .lt. work(i,7)) maxt = work(i,7) |
| if(mins .gt. work(i,6)) mins = work(i,6) |
| if(mint .gt. work(i,7)) mint = work(i,7) |
| difs = abs(work(i,6) - season(i)) |
| dift = abs(work(i,7) - trend(i)) |
| if(maxds .lt. difs) maxds = difs |
| if(maxdt .lt. dift) maxdt = dift |
| 137 continue |
| if((maxds/(maxs-mins) .lt. 0.01d0) .and. |
| & (maxdt/(maxt-mint) .lt. 0.01d0)) goto 300 |
| continue |
| j=j+1 |
| goto 100 |
| endif |
| C end Loop |
| 300 continue |
| |
| else |
| c .not. robust |
| |
| do 150 i = 1,n |
| rw(i) = 1.0d0 |
| 150 continue |
| endif |
| |
| return |
| end |
| |
| subroutine psort(a,n,ind,ni) |
| c |
| c Partial Sorting ; used for Median (MAD) computation only |
| c |
| c implicit none |
| c Arg |
| integer n,ni |
| double precision a(n) |
| integer ind(ni) |
| c Var |
| integer indu(16),indl(16),iu(16),il(16),p,jl,ju,i,j,m,k,ij,l |
| double precision t,tt |
| |
| if(n .lt. 0 .or. ni .lt. 0) return |
| |
| if(n .lt. 2 .or. ni .eq. 0) return |
| |
| jl = 1 |
| ju = ni |
| indl(1) = 1 |
| indu(1) = ni |
| i = 1 |
| j = n |
| m = 1 |
| |
| c Outer Loop |
| 161 continue |
| if(i .lt. j) go to 10 |
| |
| c _Loop_ |
| 166 continue |
| m = m-1 |
| if(m .eq. 0) return |
| i = il(m) |
| j = iu(m) |
| jl = indl(m) |
| ju = indu(m) |
| if(.not.(jl .le. ju)) goto 166 |
| |
| c while (j - i > 10) |
| 173 if(.not.(j-i .gt. 10)) goto 174 |
| |
| 10 k = i |
| ij = (i+j)/2 |
| t = a(ij) |
| if(a(i) .gt. t) then |
| a(ij) = a(i) |
| a(i) = t |
| t = a(ij) |
| endif |
| l = j |
| if(a(j) .lt. t) then |
| a(ij) = a(j) |
| a(j) = t |
| t = a(ij) |
| if(a(i) .gt. t) then |
| a(ij) = a(i) |
| a(i) = t |
| t = a(ij) |
| endif |
| endif |
| |
| 181 continue |
| l = l-1 |
| if(a(l) .le. t)then |
| tt = a(l) |
| 186 continue |
| k = k+1 |
| if(.not.(a(k) .ge. t)) goto 186 |
| |
| if(k .gt. l) goto 183 |
| |
| a(l) = a(k) |
| a(k) = tt |
| endif |
| goto 181 |
| |
| 183 continue |
| indl(m) = jl |
| indu(m) = ju |
| p = m |
| m = m+1 |
| if(l-i .le. j-k) then |
| il(p) = k |
| iu(p) = j |
| j = l |
| |
| 193 continue |
| if(jl .gt. ju) goto 166 |
| if(ind(ju) .gt. j) then |
| ju = ju-1 |
| goto 193 |
| endif |
| indl(p) = ju+1 |
| else |
| il(p) = i |
| iu(p) = l |
| i = k |
| |
| 200 continue |
| if(jl .gt. ju) goto 166 |
| if(ind(jl) .lt. i) then |
| jl = jl+1 |
| goto 200 |
| endif |
| indu(p) = jl-1 |
| endif |
| |
| goto 173 |
| c end while |
| 174 continue |
| |
| if(i .ne. 1) then |
| i = i-1 |
| 209 continue |
| i = i+1 |
| if(i .eq. j) goto 166 |
| t = a(i+1) |
| if(a(i) .gt. t) then |
| k = i |
| c repeat |
| 216 continue |
| a(k+1) = a(k) |
| k = k-1 |
| if(.not.(t .ge. a(k))) goto 216 |
| c until t >= a(k) |
| a(k+1) = t |
| endif |
| goto 209 |
| |
| endif |
| |
| goto 161 |
| c End Outer Loop |
| |
| end |
