# to unbundle, sh this file (in an empty directory) echo Makefile 1>&2 sed >Makefile <<'//GO.SYSIN DD Makefile' 's/^-//' -O= sample.o stl.o -out: a.out co2 - a.out < co2 -a.out: $O - f77 $O //GO.SYSIN DD Makefile echo co2 1>&2 sed >co2 <<'//GO.SYSIN DD co2' 's/^-//' -315.58 316.39 316.79 317.82 -318.39 318.22 316.68 315.01 -314.02 313.55 315.02 315.75 -316.52 317.1 317.79 319.22 -320.08 319.7 318.27 315.99 -314.24 314.05 315.05 316.23 -316.92 317.76 318.54 319.49 -320.64 319.85 318.7 316.96 -315.17 315.47 316.19 317.17 -318.12 318.72 319.79 320.68 -321.28 320.89 319.79 317.56 -316.46 315.59 316.85 317.87 -318.87 319.25 320.13 321.49 -322.34 321.62 319.85 317.87 -316.36 316.24 317.13 318.46 -319.57 320.23 320.89 321.54 -322.2 321.9 320.42 318.6 -316.73 317.15 317.94 318.91 -319.73 320.78 321.23 322.49 -322.59 322.35 321.61 319.24 -318.23 317.76 319.36 319.5 -320.35 321.4 322.22 323.45 -323.8 323.5 322.16 320.09 -318.26 317.66 319.47 320.7 -322.06 322.23 322.78 324.1 -324.63 323.79 322.34 320.73 -319 318.99 320.41 321.68 -322.3 322.89 323.59 324.65 -325.3 325.15 323.88 321.8 -319.99 319.86 320.88 322.36 -323.59 324.23 325.34 326.33 -327.03 326.24 325.39 323.16 -321.87 321.31 322.34 323.74 -324.61 325.58 326.55 327.81 -327.82 327.53 326.29 324.66 -323.12 323.09 324.01 325.1 -326.12 326.62 327.16 327.94 -329.15 328.79 327.53 325.65 -323.6 323.78 325.13 326.26 -326.93 327.84 327.96 329.93 -330.25 329.24 328.13 326.42 -324.97 325.29 326.56 327.73 -328.73 329.7 330.46 331.7 -332.66 332.22 331.02 329.39 -327.58 327.27 328.3 328.81 -329.44 330.89 331.62 332.85 -333.29 332.44 331.35 329.58 -327.58 327.55 328.56 329.73 -330.45 330.98 331.63 332.88 -333.63 333.53 331.9 330.08 -328.59 328.31 329.44 330.64 -331.62 332.45 333.36 334.46 -334.84 334.29 333.04 330.88 -329.23 328.83 330.18 331.5 -332.8 333.22 334.54 335.82 -336.45 335.97 334.65 332.4 -331.28 330.73 332.05 333.54 -334.65 335.06 336.32 337.39 -337.66 337.56 336.24 334.39 -332.43 332.22 333.61 334.78 -335.88 336.43 337.61 338.53 -339.06 338.92 337.39 335.72 -333.64 333.65 335.07 336.53 -337.82 338.19 339.89 340.56 -341.22 340.92 339.26 337.27 -335.66 335.54 336.71 337.79 -338.79 340.06 340.93 342.02 -342.65 341.8 340.01 337.94 -336.17 336.28 337.76 339.05 -340.18 341.04 342.16 343.01 -343.64 342.91 341.72 339.52 -337.75 337.68 339.14 340.37 -341.32 342.45 343.05 344.91 -345.77 345.3 343.98 342.41 -339.89 340.03 341.19 342.87 -343.74 344.55 345.28 347 -347.37 346.74 345.36 343.19 -340.97 341.2 342.76 343.96 -344.82 345.82 347.24 348.09 -348.66 347.9 346.27 344.21 -342.88 342.58 343.99 345.31 -345.98 346.72 347.63 349.24 -349.83 349.1 347.52 345.43 -344.48 343.89 345.29 346.54 -347.66 348.07 349.12 350.55 -351.34 350.8 349.1 347.54 -346.2 346.2 347.44 348.67 //GO.SYSIN DD co2 echo sample.f 1>&2 sed >sample.f <<'//GO.SYSIN DD sample.f' 's/^-//' -c data are 348 monthly values (29 years) in file co2 - real y(348), season(348), trend(348), rw(348), work(372,7) - logical robust - read(5,*)(y(i), i=1, 348) - n=348 - np=12 - ns=35 - nt=19 - nl=13 - no=2 - ni=1 - nsjump=4 - ntjump=2 - nljump=2 - isdeg=1 - itdeg=1 - ildeg=1 - call stl(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump, - &nljump,ni,no,rw,season, trend, work) - write(6,*) (season(i), i=1,348) - write(6,*) (trend(i), i=1,348) - write(6,*) (rw(i), i=1,348) - robust=.true. - call stlez(y,n,np,ns,isdeg,itdeg,robust,nconv,rw,season,trend, - &work) - write(6,*) nconv - write(6,*) (season(i), i=1,348) - write(6,*) (trend(i), i=1,348) - write(6,*) (rw(i), i=1,348) - end //GO.SYSIN DD sample.f echo stl.3 1>&2 sed >stl.3 <<'//GO.SYSIN DD stl.3' 's/^-//' -.SA 0 -.PH "'STL'Seasonal-Trend Decomposition'STL'" -.P -PURPOSE -.P -STL decomposes a time series into seasonal and trend components. -It returns the components and robustness weights. -.nf -.P -SYNOPSIS -.P -stl(y, n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, - nljump, ni, no, rw, season, trend, work) -integer n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, - nljump, ni, no -real y(n), rw(n), season(n), trend(n), work(n+2*np,5) -.fi -.P -ARGUMENTS -.P -.in +8 -.ti -8 -y input, time series to be decomposed. -.ti -8 -n input, number of values in y. -.ti -8 -np input, the period of the seasonal component. For example, if the -time series is monthly with a yearly cycle, then np=12. -.ti -8 -ns input, length of the seasonal smoother. -The value of ns should be an odd integer greater than or equal to 3; -ns>6 is recommended. -As ns increases -the values of the seasonal component at a given point in the seasonal cycle -(e.g., January values of a monthly series with a yearly cycle) -become smoother. -.ti -8 -nt input, length of the trend smoother. -The value of nt should be an odd integer greater than or equal to 3; -a value of nt between 1.5*np and 2*np is recommended. -As nt increases the values of the trend component become smoother. -.ti -8 -nl input, length of the low-pass filter. -The value of nl should be an odd integer greater than or equal to 3; -the smallest odd integer greater than or equal to np is recommended. -.ti -8 -isdeg input, degree of locally-fitted polynomial in seasonal smoothing. -The value is 0 or 1. -.ti -8 -itdeg input, degree of locally-fitted polynomial in trend smoothing. -The value is 0 or 1. -.ti -8 -ildeg input, degree of locally-fitted polynomial in low-pass smoothing. -The value is 0 or 1. -.ti -8 -nsjump input, skipping value for seasonal smoothing. -The seasonal smoother skips ahead nsjump points and then linearly interpolates in between. -The value of nsjump should be a positive integer; -if nsjump=1, a seasonal smooth is calculated at all n points. -To make the procedure run faster, -a reasonable choice for nsjump is 10%-20% of ns. -.ti -8 -ntjump input, skipping value for trend smoothing. -.ti -8 -nljump input, skipping value for the low-pass filter. -.ti -8 -ni input, number of loops for updating the seasonal and trend components. -The value of ni should be a positive integer. -See the next argument for advice on the choice of ni. -.ti -8 -no input, number of iterations of robust fitting. -The value of no should be a nonnegative integer. -If the data are well behaved without outliers, -then robustness iterations are not needed. -In this case set no=0, and set ni=2 to 5 -depending on how much security you want that the seasonal-trend looping converges. -If outliers are present then no=3 is a very secure value -unless the outliers are radical, -in which case no=5 or even 10 might be better. -If no>0 then set ni to 1 or 2. -.ti -8 -rw output, final robustness weights. All rw are 1 if no=0. -.ti -8 -season output, seasonal component. -.ti -8 -trend output, trend component. -.ti -8 -work workspace of (n+2*np)*5 locations. -.in -8 //GO.SYSIN DD stl.3 echo stl.doc 1>&2 sed >stl.doc <<'//GO.SYSIN DD stl.doc' 's/^-//' - PURPOSE - STL decomposes a time series into seasonal and trend components. - It returns the components and robustness weights. - SYNOPSIS - stl(y, n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, - nljump, ni, no, rw, season, trend, work) - integer n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, - nljump, ni, no - real y(n), rw(n), season(n), trend(n), work(n+2*np,5) - ARGUMENTS - y input, time series to be decomposed. - n input, number of values in y. - np input, the period of the seasonal component. For example, - if the time series is monthly with a yearly cycle, then - np=12. - ns input, length of the seasonal smoother. The value of ns - should be an odd integer greater than or equal to 3; ns>6 - is recommended. As ns increases the values of the - seasonal component at a given point in the seasonal cycle - (e.g., January values of a monthly series with a yearly - cycle) become smoother. - nt input, length of the trend smoother. The value of nt - should be an odd integer greater than or equal to 3; a - value of nt between 1.5*np and 2*np is recommended. As - nt increases the values of the trend component become - smoother. - nl input, length of the low-pass filter. The value of nl - should be an odd integer greater than or equal to 3; the - smallest odd integer greater than or equal to np is - recommended. - isdeg input, degree of locally-fitted polynomial in seasonal - smoothing. The value is 0 or 1. - itdeg input, degree of locally-fitted polynomial in trend - smoothing. The value is 0 or 1. - ildeg input, degree of locally-fitted polynomial in low-pass - smoothing. The value is 0 or 1. - nsjump input, skipping value for seasonal smoothing. The - seasonal smoother skips ahead nsjump points and then - linearly interpolates in between. The value of nsjump - should be a positive integer; if nsjump=1, a seasonal - smooth is calculated at all n points. To make the - procedure run faster, a reasonable choice for nsjump is - 10%-20% of ns. - ntjump input, skipping value for trend smoothing. - nljump input, skipping value for the low-pass filter. - ni input, number of loops for updating the seasonal and - trend components. The value of ni should be a positive - integer. See the next argument for advice on the choice - of ni. - no input, number of iterations of robust fitting. The value - of no should be a nonnegative integer. If the data are - well behaved without outliers, then robustness iterations - are not needed. In this case set no=0, and set ni=2 to 5 - depending on how much security you want that the - seasonal-trend looping converges. If outliers are - present then no=3 is a very secure value unless the - outliers are radical, in which case no=5 or even 10 might - be better. If no>0 then set ni to 1 or 2. - rw output, final robustness weights. All rw are 1 if no=0. - season output, seasonal component. - trend output, trend component. - work workspace of (n+2*np)*5 locations. //GO.SYSIN DD stl.doc echo stl.f 1>&2 sed >stl.f <<'//GO.SYSIN DD stl.f' 's/^-//' - subroutine stl(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump, -& nljump,ni,no,rw,season,trend,work) - integer n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, -& nljump, ni, no, k - integer newns, newnt, newnl, newnp - real y(n), rw(n), season(n), trend(n), work(n+2*np,5) - logical userw - userw = .false. - k = 0 - do 23000 i = 1,n - trend(i) = 0.0 -23000 continue - newns = max0(3,ns) - newnt = max0(3,nt) - newnl = max0(3,nl) - newnp = max0(2,np) - if(.not.(mod(newns,2) .eq. 0))goto 23002 - newns = newns + 1 -23002 continue - if(.not.(mod(newnt,2) .eq. 0))goto 23004 - newnt = newnt + 1 -23004 continue - if(.not.(mod(newnl,2) .eq. 0))goto 23006 - newnl = newnl + 1 -23006 continue -23008 continue - call onestp(y,n,newnp,newns,newnt,newnl,isdeg,itdeg,ildeg,nsjump, -& ntjump,nljump,ni,userw,rw,season, trend, work) - k = k+1 - if(.not.(k .gt. no))goto 23011 - goto 23010 -23011 continue - do 23013 i = 1,n - work(i,1) = trend(i)+season(i) -23013 continue - call rwts(y,n,work(1,1),rw) - userw = .true. -23009 goto 23008 -23010 continue - if(.not.(no .le. 0))goto 23015 - do 23017 i = 1,n - rw(i) = 1.0 -23017 continue -23015 continue - return - end - subroutine ess(y,n,len,ideg,njump,userw,rw,ys,res) - integer n, len, ideg, njump, newnj, nleft, nright, nsh, k, i, j - real y(n), rw(n), ys(n), res(n), delta - logical ok, userw - if(.not.(n .lt. 2))goto 23019 - ys(1) = y(1) - return -23019 continue - newnj = min0(njump, n-1) - if(.not.(len .ge. n))goto 23021 - nleft = 1 - nright = n - do 23023 i = 1,n,newnj - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(.not.( .not. ok))goto 23025 - ys(i) = y(i) -23025 continue -23023 continue - goto 23022 -23021 continue - if(.not.(newnj .eq. 1))goto 23027 - nsh = (len+1)/2 - nleft = 1 - nright = len - do 23029 i = 1,n - if(.not.(i .gt. nsh .and. nright .ne. n))goto 23031 - nleft = nleft+1 - nright = nright+1 -23031 continue - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(.not.( .not. ok))goto 23033 - ys(i) = y(i) -23033 continue -23029 continue - goto 23028 -23027 continue - nsh = (len+1)/2 - do 23035 i = 1,n,newnj - if(.not.(i .lt. nsh))goto 23037 - nleft = 1 - nright = len - goto 23038 -23037 continue - if(.not.(i .ge. n-nsh+1))goto 23039 - nleft = n-len+1 - nright = n - goto 23040 -23039 continue - nleft = i-nsh+1 - nright = len+i-nsh -23040 continue -23038 continue - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(.not.( .not. ok))goto 23041 - ys(i) = y(i) -23041 continue -23035 continue -23028 continue -23022 continue - if(.not.(newnj .ne. 1))goto 23043 - do 23045 i = 1,n-newnj,newnj - delta = (ys(i+newnj)-ys(i))/float(newnj) - do 23047 j = i+1,i+newnj-1 - ys(j) = ys(i)+delta*float(j-i) -23047 continue -23045 continue - k = ((n-1)/newnj)*newnj+1 - if(.not.(k .ne. n))goto 23049 - call est(y,n,len,ideg,float(n),ys(n),nleft,nright,res,userw,rw,ok) - if(.not.( .not. ok))goto 23051 - ys(n) = y(n) -23051 continue - if(.not.(k .ne. n-1))goto 23053 - delta = (ys(n)-ys(k))/float(n-k) - do 23055 j = k+1,n-1 - ys(j) = ys(k)+delta*float(j-k) -23055 continue -23053 continue -23049 continue -23043 continue - return - end - subroutine est(y,n,len,ideg,xs,ys,nleft,nright,w,userw,rw,ok) - integer n, len, ideg, nleft, nright, j - real y(n), w(n), rw(n), xs, ys, range, h, h1, h9, a, b, c, r - logical userw,ok - range = float(n)-float(1) - h = amax1(xs-float(nleft),float(nright)-xs) - if(.not.(len .gt. n))goto 23057 - h = h+float((len-n)/2) -23057 continue - h9 = .999*h - h1 = .001*h - a = 0.0 - do 23059 j = nleft,nright - w(j) = 0. - r = abs(float(j)-xs) - if(.not.(r .le. h9))goto 23061 - if(.not.(r .le. h1))goto 23063 - w(j) = 1. - goto 23064 -23063 continue - w(j) = (1.0-(r/h)**3)**3 -23064 continue - if(.not.(userw))goto 23065 - w(j) = rw(j)*w(j) -23065 continue - a = a+w(j) -23061 continue -23059 continue - if(.not.(a .le. 0.0))goto 23067 - ok = .false. - goto 23068 -23067 continue - ok = .true. - do 23069 j = nleft,nright - w(j) = w(j)/a -23069 continue - if(.not.((h .gt. 0.) .and. (ideg .gt. 0)))goto 23071 - a = 0.0 - do 23073 j = nleft,nright - a = a+w(j)*float(j) -23073 continue - b = xs-a - c = 0.0 - do 23075 j = nleft,nright - c = c+w(j)*(float(j)-a)**2 -23075 continue - if(.not.(sqrt(c) .gt. .001*range))goto 23077 - b = b/c - do 23079 j = nleft,nright - w(j) = w(j)*(b*(float(j)-a)+1.0) -23079 continue -23077 continue -23071 continue - ys = 0.0 - do 23081 j = nleft,nright - ys = ys+w(j)*y(j) -23081 continue -23068 continue - return - end - subroutine fts(x,n,np,trend,work) - integer n, np - real x(n), trend(n), work(n) - call ma(x,n,np,trend) - call ma(trend,n-np+1,np,work) - call ma(work,n-2*np+2,3,trend) - return - end - subroutine ma(x, n, len, ave) - integer n, len, i, j, k, m, newn - real x(n), ave(n), flen, v - newn = n-len+1 - flen = float(len) - v = 0.0 - do 23083 i = 1,len - v = v+x(i) -23083 continue - ave(1) = v/flen - if(.not.(newn .gt. 1))goto 23085 - k = len - m = 0 - do 23087 j = 2, newn - k = k+1 - m = m+1 - v = v-x(m)+x(k) - ave(j) = v/flen -23087 continue -23085 continue - return - end - subroutine onestp(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump, -& nljump,ni,userw,rw,season,trend,work) - integer n,ni,np,ns,nt,nsjump,ntjump,nl,nljump,isdeg,itdeg,ildeg - real y(n),rw(n),season(n),trend(n),work(n+2*np,5) - logical userw - do 23089 j = 1,ni - do 23091 i = 1,n - work(i,1) = y(i)-trend(i) -23091 continue - call ss(work(1,1),n,np,ns,isdeg,nsjump,userw,rw,work(1,2),work(1, -& 3),work(1,4),work(1,5),season) - call fts(work(1,2),n+2*np,np,work(1,3),work(1,1)) - call ess(work(1,3),n,nl,ildeg,nljump,.false.,work(1,4),work(1,1), -& work(1,5)) - do 23093 i = 1,n - season(i) = work(np+i,2)-work(i,1) -23093 continue - do 23095 i = 1,n - work(i,1) = y(i)-season(i) -23095 continue - call ess(work(1,1),n,nt,itdeg,ntjump,userw,rw,trend,work(1,3)) -23089 continue - return - end - subroutine rwts(y,n,fit,rw) - integer mid(2), n - real y(n), fit(n), rw(n), cmad, c9, c1, r - do 23097 i = 1,n - rw(i) = abs(y(i)-fit(i)) -23097 continue - mid(1) = n/2+1 - mid(2) = n-mid(1)+1 - call psort(rw,n,mid,2) - cmad = 3.0*(rw(mid(1))+rw(mid(2))) - c9 = .999*cmad - c1 = .001*cmad - do 23099 i = 1,n - r = abs(y(i)-fit(i)) - if(.not.(r .le. c1))goto 23101 - rw(i) = 1. - goto 23102 -23101 continue - if(.not.(r .le. c9))goto 23103 - rw(i) = (1.0-(r/cmad)**2)**2 - goto 23104 -23103 continue - rw(i) = 0. -23104 continue -23102 continue -23099 continue - return - end - subroutine ss(y,n,np,ns,isdeg,nsjump,userw,rw,season,work1,work2, -& work3,work4) - integer n, np, ns, isdeg, nsjump, nright, nleft, i, j, k - real y(n), rw(n), season(n+2*np), work1(n), work2(n), work3(n), -& work4(n), xs - logical userw,ok - j=1 -23105 if(.not.(j .le. np))goto 23107 - k = (n-j)/np+1 - do 23108 i = 1,k - work1(i) = y((i-1)*np+j) -23108 continue - if(.not.(userw))goto 23110 - do 23112 i = 1,k - work3(i) = rw((i-1)*np+j) -23112 continue -23110 continue - call ess(work1,k,ns,isdeg,nsjump,userw,work3,work2(2),work4) - xs = 0 - nright = min0(ns,k) - call est(work1,k,ns,isdeg,xs,work2(1),1,nright,work4,userw,work3, -& ok) - if(.not.( .not. ok))goto 23114 - work2(1) = work2(2) -23114 continue - xs = k+1 - nleft = max0(1,k-ns+1) - call est(work1,k,ns,isdeg,xs,work2(k+2),nleft,k,work4,userw,work3, -& ok) - if(.not.( .not. ok))goto 23116 - work2(k+2) = work2(k+1) -23116 continue - do 23118 m = 1,k+2 - season((m-1)*np+j) = work2(m) -23118 continue - j=j+1 - goto 23105 -23107 continue - return - end - subroutine stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, -& season, trend, work) - logical robust - integer n, i, j, np, ns, no, nt, nl, ni, nsjump, ntjump, nljump, -& newns, newnp - integer isdeg, itdeg, ildeg - real y(n), rw(n), season(n), trend(n), work(n+2*np,7) - real maxs, mins, maxt, mint, maxds, maxdt, difs, dift - ildeg = itdeg - newns = max0(3,ns) - if(.not.(mod(newns,2) .eq. 0))goto 23120 - newns = newns+1 -23120 continue - newnp = max0(2,np) - nt = (1.5*newnp)/(1 - 1.5/newns) + 0.5 - nt = max0(3,nt) - if(.not.(mod(nt,2) .eq. 0))goto 23122 - nt = nt+1 -23122 continue - nl = newnp - if(.not.(mod(nl,2) .eq. 0))goto 23124 - nl = nl+1 -23124 continue - if(.not.(robust))goto 23126 - ni = 1 - goto 23127 -23126 continue - ni = 2 -23127 continue - 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 23128 i = 1,n - trend(i) = 0.0 -23128 continue - call onestp(y,n,newnp,newns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump, -& nljump,ni,.false.,rw,season,trend,work) - no = 0 - if(.not.(robust))goto 23130 - j=1 -23132 if(.not.(j .le. 15))goto 23134 - do 23135 i = 1,n - work(i,6) = season(i) - work(i,7) = trend(i) - work(i,1) = trend(i)+season(i) -23135 continue - call rwts(y,n,work(1,1),rw) - call onestp(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 23137 i = 2,n - if(.not.(maxs .lt. work(i,6)))goto 23139 - maxs = work(i,6) -23139 continue - if(.not.(maxt .lt. work(i,7)))goto 23141 - maxt = work(i,7) -23141 continue - if(.not.(mins .gt. work(i,6)))goto 23143 - mins = work(i,6) -23143 continue - if(.not.(mint .gt. work(i,7)))goto 23145 - mint = work(i,7) -23145 continue - difs = abs(work(i,6) - season(i)) - dift = abs(work(i,7) - trend(i)) - if(.not.(maxds .lt. difs))goto 23147 - maxds = difs -23147 continue - if(.not.(maxdt .lt. dift))goto 23149 - maxdt = dift -23149 continue -23137 continue - if(.not.((maxds/(maxs-mins) .lt. .01) .and. (maxdt/(maxt-mint) -& .lt. .01)))goto 23151 - goto 23134 -23151 continue - j=j+1 - goto 23132 -23134 continue -23130 continue - if(.not.( .not. robust))goto 23153 - do 23155 i = 1,n - rw(i) = 1.0 -23155 continue -23153 continue - return - end - subroutine psort(a,n,ind,ni) - real a(n) - integer n,ind(ni),ni - integer indu(16),indl(16),iu(16),il(16),p,jl,ju,i,j,m,k,ij,l - real t,tt - if(.not.(n .lt. 0 .or. ni .lt. 0))goto 23157 - return -23157 continue - if(.not.(n .lt. 2 .or. ni .eq. 0))goto 23159 - return -23159 continue - jl = 1 - ju = ni - indl(1) = 1 - indu(1) = ni - i = 1 - j = n - m = 1 -23161 continue - if(.not.(i .lt. j))goto 23164 - go to 10 -23164 continue -23166 continue - m = m-1 - if(.not.(m .eq. 0))goto 23169 - goto 23163 -23169 continue - i = il(m) - j = iu(m) - jl = indl(m) - ju = indu(m) - if(.not.(jl .le. ju))goto 23171 -23173 if(.not.(j-i .gt. 10))goto 23174 -10 k = i - ij = (i+j)/2 - t = a(ij) - if(.not.(a(i) .gt. t))goto 23175 - a(ij) = a(i) - a(i) = t - t = a(ij) -23175 continue - l = j - if(.not.(a(j) .lt. t))goto 23177 - a(ij) = a(j) - a(j) = t - t = a(ij) - if(.not.(a(i) .gt. t))goto 23179 - a(ij) = a(i) - a(i) = t - t = a(ij) -23179 continue -23177 continue -23181 continue - l = l-1 - if(.not.(a(l) .le. t))goto 23184 - tt = a(l) -23186 continue - k = k+1 -23187 if(.not.(a(k) .ge. t))goto 23186 - if(.not.(k .gt. l))goto 23189 - goto 23183 -23189 continue - a(l) = a(k) - a(k) = tt -23184 continue -23182 goto 23181 -23183 continue - indl(m) = jl - indu(m) = ju - p = m - m = m+1 - if(.not.(l-i .le. j-k))goto 23191 - il(p) = k - iu(p) = j - j = l -23193 continue - if(.not.(jl .gt. ju))goto 23196 - goto 23167 -23196 continue - if(.not.(ind(ju) .le. j))goto 23198 - goto 23195 -23198 continue - ju = ju-1 -23194 goto 23193 -23195 continue - indl(p) = ju+1 - goto 23192 -23191 continue - il(p) = i - iu(p) = l - i = k -23200 continue - if(.not.(jl .gt. ju))goto 23203 - goto 23167 -23203 continue - if(.not.(ind(jl) .ge. i))goto 23205 - goto 23202 -23205 continue - jl = jl+1 -23201 goto 23200 -23202 continue - indu(p) = jl-1 -23192 continue - goto 23173 -23174 continue - if(.not.(i .eq. 1))goto 23207 - goto 23168 -23207 continue - i = i-1 -23209 continue - i = i+1 - if(.not.(i .eq. j))goto 23212 - goto 23211 -23212 continue - t = a(i+1) - if(.not.(a(i) .gt. t))goto 23214 - k = i -23216 continue - a(k+1) = a(k) - k = k-1 -23217 if(.not.(t .ge. a(k)))goto 23216 - a(k+1) = t -23214 continue -23210 goto 23209 -23211 continue -23171 continue -23167 goto 23166 -23168 continue -23162 goto 23161 -23163 continue - return - end //GO.SYSIN DD stl.f echo stl.int 1>&2 sed >stl.int <<'//GO.SYSIN DD stl.int' 's/^-//' - Routines for the STL Seasonal-Trend Decomposition Procedure - - - - This directory contains routines for carrying out the STL - decomposition procedure on seasonal time series. Most users - will need only to call one or both of the top-level routines - stl and stlez. They are documented in files stl.doc and - stlez.doc. Routine stl allows specification of all - parameters; routine stlez (pronounced "S-T-L easy") selects - default values for most parameters and carries out - robustness iterations to convergence. The STL statistical - method can be applied in principle even when there are - missing values in the time series, but these routines do not - allow missing values. - - On a UNIX system, just type "make" to compile *.f and run - the sample program. On a PC, change the names of the - Fortran files by adding "or" to each, compile *.for, and - load; files renam,bat, compile.bat, load.bat, and objects - may be helpful. - - Estimation of Other Components - - In some applications components other than the seasonal and - trend need to be estimated. Trading-day estimation for - economic time series is one example. The routines supplied - here have been constructed so that Fortran or Ratfor code - that carries out estimation of such components can be - inserted easily. - - An illustration is given by the following modification of - the Ratfor routine in stl.r. In this modification a routine - called "yours" estimates a component called "other". Notice - that the robustness weights enter yours; that is, if stl is - carrying out robust estimation of the seasonal and trend - components then it is appropriate to robustly estimate other - components. Usually the routine yours would need workspace, - though none is given here. - - subroutine stl(y, n, np, ns, nt, nl, isdeg, itdeg, ildeg, - nsjump, ntjump, nljump, ni, no, rw, - season, trend, other, more, work) - integer n, np, ns, nt, nl, nsjump, ntjump, nljump, ni, no, k - integer newns, newnt, newnl, newnp - real y(n), rw(n), season(n), trend(n), other(n), more(n), - work(n+2*np,5) - logical userw - userw = .false. - k = 0 - do i = 1, n { - trend(i) = 0.0 - other(i) = 0.0 - } - newns = max0(3,ns) - newnt = max0(3,nt) - newnl = max0(3,nl) - newnp = max0(2,np) - if(mod(newns,2)==0) newns = newns + 1 #make odd - if(mod(newnt,2)==0) newnt = newnt + 1 - if(mod(newnl,2)==0) newnl = newnl + 1 - repeat { - do i = 1, n - more(i) = y(i) - other(i) - call onestp(more, n, newnp, newns, newnt, newnl, isdeg, - itdeg, ildeg, nsjump, ntjump, nljump, ni, - userw, rw, season, trend, work) - do i = 1, n - more(i) = y(i) - season(i) - trend(i) - call yours(more, n, userw, rw, other) - k = k+1 - if(k>no) break - do i = 1, n - work(i,1) = trend(i)+season(i)+other(i) - call rwts(y, n, work(1, 1), rw) - userw = .true. - } - if (no<=0) - do i = 1, n - rw(i) = 1.0 - return - end - - - - - - - - - - - - - - - - - - - - //GO.SYSIN DD stl.int echo stl.r 1>&2 sed >stl.r <<'//GO.SYSIN DD stl.r' 's/^-//' -subroutine stl(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,ni,no,rw, - season,trend,work) - -integer n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ni, no, k -integer newns, newnt, newnl, newnp -real y(n), rw(n), season(n), trend(n), work(n+2*np,5) -logical userw - -userw = .false. -k = 0 -do i = 1,n - trend(i) = 0.0 -newns = max0(3,ns) -newnt = max0(3,nt) -newnl = max0(3,nl) -newnp = max0(2,np) -if(mod(newns,2)==0) newns = newns + 1 #make odd -if(mod(newnt,2)==0) newnt = newnt + 1 -if(mod(newnl,2)==0) newnl = newnl + 1 -repeat { - call onestp(y,n,newnp,newns,newnt,newnl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,n i,userw,rw,season, trend, work) - k = k+1 - if(k>no) break - do i = 1,n - work(i,1) = trend(i)+season(i) - call rwts(y,n,work(1,1),rw) - userw = .true. - } -if (no<=0) - do i = 1,n - rw(i) = 1.0 -return -end - - - -subroutine ess(y,n,len,ideg,njump,userw,rw,ys,res) - -integer n, len, ideg, njump, newnj, nleft, nright, nsh, k, i, j -real y(n), rw(n), ys(n), res(n), delta -logical ok, userw - -if(n<2) { ys(1) = y(1); return } -newnj = min0(njump, n-1) -if(len>=n) { #len > or = n - nleft = 1 - nright = n - do i = 1,n,newnj { - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(!ok) ys(i) = y(i) - } - } -else if(newnj==1) { # newnj equal to one, len less than n - nsh = (len+1)/2 - nleft = 1 - nright = len - do i = 1,n { # fitted value at i - if(i>nsh && nright!=n) { - nleft = nleft+1 - nright = nright+1 - } - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(!ok) ys(i) = y(i) - } - } -else { # newnj greater than one, len less than n - nsh = (len+1)/2 - do i = 1,n,newnj { # fitted value at i - if(i=n-nsh+1) { - nleft = n-len+1 - nright = n - } - else { - nleft = i-nsh+1 - nright = len+i-nsh - } - call est(y,n,len,ideg,float(i),ys(i),nleft,nright,res,userw,rw,ok) - if(!ok) ys(i) = y(i) - } - } -if(newnj!=1){ - do i = 1,n-newnj,newnj { - delta = (ys(i+newnj)-ys(i))/float(newnj) - do j = i+1,i+newnj-1 - ys(j) = ys(i)+delta*float(j-i) - } - k = ((n-1)/newnj)*newnj+1 - if(k!=n) { - call est(y,n,len,ideg,float(n),ys(n),nleft,nright,res,userw,rw,ok) - if(!ok) ys(n) = y(n) - if(k!=n-1) { - delta = (ys(n)-ys(k))/float(n-k) - do j = k+1,n-1 - ys(j) = ys(k)+delta*float(j-k) - } - } - } -return -end - - -subroutine est(y,n,len,ideg,xs,ys,nleft,nright,w,userw,rw,ok) - -integer n, len, ideg, nleft, nright, j -real y(n), w(n), rw(n), xs, ys, range, h, h1, h9, a, b, c, r -logical userw,ok - -range = float(n)-float(1) -h = amax1(xs-float(nleft),float(nright)-xs) -if (len>n) h = h+float((len-n)/2) -h9 = .999*h -h1 = .001*h -# compute weights -a = 0.0 -do j = nleft,nright { - w(j) = 0. - r = abs(float(j)-xs) - if (r<=h9) { - if (r<=h1) w(j) = 1. - else w(j) = (1.0-(r/h)**3)**3 - if (userw) w(j) = rw(j)*w(j) - a = a+w(j) - } - } -if (a<=0.0) - ok = .false. -else { # weighted least squares - ok = .true. - do j = nleft,nright # make sum of w(j) == 1 - w(j) = w(j)/a - if ((h>0.)&(ideg>0)) { # use linear fit - a = 0.0 - do j = nleft,nright # weighted center of x values - a = a+w(j)*float(j) - b = xs-a - c = 0.0 - do j = nleft,nright - c = c+w(j)*(float(j)-a)**2 - if (sqrt(c)>.001*range) { - b = b/c -# points are spread out enough to compute slope - do j = nleft,nright - w(j) = w(j)*(b*(float(j)-a)+1.0) - } - } - ys = 0.0 - do j = nleft,nright - ys = ys+w(j)*y(j) - } -return -end - - - -subroutine fts(x,n,np,trend,work) - -integer n, np -real x(n), trend(n), work(n) - -call ma(x,n,np,trend) -call ma(trend,n-np+1,np,work) -call ma(work,n-2*np+2,3,trend) -return -end - - -subroutine ma(x, n, len, ave) - -integer n, len, i, j, k, m, newn -real x(n), ave(n), flen, v - -newn = n-len+1 -flen = float(len) -v = 0.0 -# get the first average -do i = 1,len - v = v+x(i) -ave(1) = v/flen -if (newn>1) { - k = len - m = 0 - do j = 2, newn { -# window down the array - k = k+1 - m = m+1 - v = v-x(m)+x(k) - ave(j) = v/flen - } -} -return -end - - -subroutine onestp(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,ni, - userw,rw,season,trend,work) - -integer n,ni,np,ns,nt,nsjump,ntjump,nl,nljump,isdeg,itdeg,ildeg -real y(n),rw(n),season(n),trend(n),work(n+2*np,5) -logical userw - -do j = 1,ni { - do i = 1,n - work(i,1) = y(i)-trend(i) - call ss(work(1,1),n,np,ns,isdeg,nsjump,userw,rw,work(1,2),work(1,3),work(1,4), - work(1,5),season) - call fts(work(1,2),n+2*np,np,work(1,3),work(1,1)) - call ess(work(1,3),n,nl,ildeg,nljump,.false.,work(1,4),work(1,1),work(1,5)) - do i = 1,n - season(i) = work(np+i,2)-work(i,1) - do i = 1,n - work(i,1) = y(i)-season(i) - call ess(work(1,1),n,nt,itdeg,ntjump,userw,rw,trend,work(1,3)) - } -return -end - - -subroutine rwts(y,n,fit,rw) - -integer mid(2), n -real y(n), fit(n), rw(n), cmad, c9, c1, r - -do i = 1,n - rw(i) = abs(y(i)-fit(i)) -mid(1) = n/2+1 -mid(2) = n-mid(1)+1 -call psort(rw,n,mid,2) -cmad = 3.0*(rw(mid(1))+rw(mid(2))) #6 * median abs resid -c9 = .999*cmad -c1 = .001*cmad -do i = 1,n { - r = abs(y(i)-fit(i)) - if (r<=c1) rw(i) = 1. - else if (r<=c9) rw(i) = (1.0-(r/cmad)**2)**2 - else rw(i) = 0. - } -return -end - - - -subroutine ss(y,n,np,ns,isdeg,nsjump,userw,rw,season,work1,work2,work3,work4) - -integer n, np, ns, isdeg, nsjump, nright, nleft, i, j, k -real y(n), rw(n), season(n+2*np), work1(n), work2(n), work3(n), work4(n), xs -logical userw,ok - -for(j=1; j<=np; j=j+1){ - k = (n-j)/np+1 - do i = 1,k - work1(i) = y((i-1)*np+j) - if(userw) - do i = 1,k - work3(i) = rw((i-1)*np+j) - call ess(work1,k,ns,isdeg,nsjump,userw,work3,work2(2),work4) - xs = 0 - nright = min0(ns,k) - call est(work1,k,ns,isdeg,xs,work2(1),1,nright,work4,userw,work3,ok) - if(!ok) work2(1) = work2(2) - xs = k+1 - nleft = max0(1,k-ns+1) - call est(work1,k,ns,isdeg,xs,work2(k+2),nleft,k,work4,userw,work3,ok) - if(!ok) work2(k+2) = work2(k+1) - do m = 1,k+2 - season((m-1)*np+j) = work2(m) - } -return -end -subroutine stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, season, trend, work) - -logical robust -integer n, i, j, np, ns, no, nt, nl, ni, nsjump, ntjump, nljump, newns, newnp -integer isdeg, itdeg, ildeg -real y(n), rw(n), season(n), trend(n), work(n+2*np,7) -real maxs, mins, maxt, mint, maxds, maxdt, difs, dift - -ildeg = itdeg -newns = max0(3,ns) -if(mod(newns,2)==0) newns = newns+1 -newnp = max0(2,np) -nt = (1.5*newnp)/(1 - 1.5/newns) + 0.5 -nt = max0(3,nt) -if(mod(nt,2)==0) nt = nt+1 -nl = newnp -if(mod(nl,2)==0) nl = nl+1 -if(robust) ni = 1 -else ni = 2 -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 i = 1,n - trend(i) = 0.0 -call onestp(y,n,newnp,newns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,nljump,ni, - .false.,rw,season,trend,work) -no = 0 -if(robust){ - for(j=1; j<=15; j=j+1){ #robustness iterations - do i = 1,n{ #initialize for testing - work(i,6) = season(i) - work(i,7) = trend(i) - work(i,1) = trend(i)+season(i) - } - call rwts(y,n,work(1,1),rw) - call onestp(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 i = 2,n{ - if(maxswork(i,6)) mins = work(i,6) - if(mint>work(i,7)) mint = work(i,7) - difs = abs(work(i,6) - season(i)) - dift = abs(work(i,7) - trend(i)) - if(maxds10) { -# first order a(i),a(j),a((i+j)/2), and use median to split the data - 10 k = i - ij = (i+j)/2 - t = a(ij) - if (a(i)>t) { - a(ij) = a(i) - a(i) = t - t = a(ij) - } - l = j - if (a(j)t) { - a(ij) = a(i) - a(i) = t - t = a(ij) - } - } - repeat { - l = l-1 - if (a(l)<=t) { - tt = a(l) - repeat -# split the data into a(i to l).lt.t, a(k to j).gt.t - k = k+1 - until(a(k)>=t) - if (k>l) break - a(l) = a(k) - a(k) = tt - } - } - indl(m) = jl - indu(m) = ju - p = m - m = m+1 -# split the larger of the segments - if (l-i<=j-k) { - il(p) = k - iu(p) = j - j = l - repeat { - if (jl>ju) next 3 - if (ind(ju)<=j) break - ju = ju-1 - } - indl(p) = ju+1 - } - else { - il(p) = i - iu(p) = l - i = k - repeat { -# skip all segments not corresponding to an entry in ind - if (jl>ju) next 3 - if (ind(jl)>=i) break - jl = jl+1 - } - indu(p) = jl-1 - } - } - if (i==1) break - i = i-1 - repeat { - i = i+1 - if (i==j) break - t = a(i+1) - if (a(i)>t) { - k = i - repeat { - a(k+1) = a(k) - k = k-1 - } - until(t>=a(k)) - a(k+1) = t - } - } - } - } - } -return -end //GO.SYSIN DD stl.r echo stlez.3 1>&2 sed >stlez.3 <<'//GO.SYSIN DD stlez.3' 's/^-//' -.SA 0 -.PH "'STLEZ'Seasonal-Trend Decomposition with Defaults'STLEZ'" -.P -PURPOSE -.P -STLEZ offers easy use of the STL procedure by defaulting most parameters values. -.nf -.P -SYNOPSIS -.P -stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, season, - trend, work) -integer n, np, ns, no isdeg, itdeg -real y(n), rw(n), season(n), trend(n), work(n+2*np,7) -logical robust -.fi -.P -ARGUMENTS -.P -.in +8 -.ti -8 -y input, time series to be decomposed. -.ti -8 -n input, number of values in y. -.ti -8 -np input, the period of the seasonal component. For example, if the -time series is monthly with a yearly cycle, then np=12. -.ti -8 -ns input, length of the seasonal smoother. -The value of ns should be an odd integer greater than or equal to 3; -ns>6 is recommended. -As ns increases -the values of the seasonal component at a given point in the seasonal cycle -(e.g., January values of a monthly series with a yearly cycle) -become smoother. -.ti -8 -isdeg input, degree of locally-fitted polynomial in seasonal smoothing. -The value is 0 or 1. -.ti -8 -itdeg input, degree of locally-fitted polynomial in trend smoothing. -The value is 0 or 1. -.ti -8 -robust input, .true. if robustness iterations are to be used, .false. otherwise. -.ti -8 -rw output, final robustness weights. -All rw are 1 if robust=.false. -.ti -8 -season output, seasonal component. -.ti -8 -trend output, trend component. -.ti -8 -no output, number of robustness iterations. The iterations end -if a convergence criterion is met or if the number is 15. -.ti -8 -work workspace of (n+2*np)*7 locations. -.in -8 -.P -REFERENCES -.P -R.B. Cleveland, W.S.Cleveland, J.E. McRae, and I. Terpenning, -STL: A Seasonal-Trend Decomposition Procedure Based on Loess, Statistics Research -Report, AT&T Bell Laboratories. -.P -STLEZ sets parameters to the values recommended in the above reference. -Thus if you call STLEZ, you will be using the values summarized below. -See documentation for STL for definitions of the parameters. -.br -nt = smallest odd integer greater than or equal to -.ce -(1.5*np) / (1-(1.5/ns)). -.br -nl = smallest odd integer greater than or equal to np. -.br -nsjump = ns/10; ntjump = nt/10; nljump = nl/10. -.br -ni = 2 if robust=.false.; else ni=1. -.br -no = 0 if robust=.false.; else robustness iterations are carried out -until convergence of both seasonal and trend components, -with 15 iterations maximum. -Convergence occurs if the maximum changes in individual seasonal and trend fits -are less than 1% of the component's range after the previous iteration. //GO.SYSIN DD stlez.3 echo stlez.doc 1>&2 sed >stlez.doc <<'//GO.SYSIN DD stlez.doc' 's/^-//' - PURPOSE - STLEZ offers easy use of the STL procedure by defaulting most - parameters values. - SYNOPSIS - stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, season, - trend, work) - integer n, np, ns, no isdeg, itdeg - real y(n), rw(n), season(n), trend(n), work(n+2*np,7) - logical robust - ARGUMENTS - y input, time series to be decomposed. - n input, number of values in y. - np input, the period of the seasonal component. For example, - if the time series is monthly with a yearly cycle, then - np=12. - ns input, length of the seasonal smoother. The value of ns - should be an odd integer greater than or equal to 3; ns>6 - is recommended. As ns increases the values of the - seasonal component at a given point in the seasonal cycle - (e.g., January values of a monthly series with a yearly - cycle) become smoother. - isdeg input, degree of locally-fitted polynomial in seasonal - smoothing. The value is 0 or 1. - itdeg input, degree of locally-fitted polynomial in trend - smoothing. The value is 0 or 1. - robust input, .true. if robustness iterations are to be used, - .false. otherwise. - rw output, final robustness weights. All rw are 1 if - robust=.false. - season output, seasonal component. - trend output, trend component. - no output, number of robustness iterations. The iterations - end if a convergence criterion is met or if the number is - 15. - work workspace of (n+2*np)*7 locations. - REFERENCES - R.B. Cleveland, W.S.Cleveland, J.E. McRae, and I. Terpenning, - STL: A Seasonal-Trend Decomposition Procedure Based on Loess, - Statistics Research Report, AT&T Bell Laboratories. - STLEZ sets parameters to the values recommended in the above - reference. Thus if you call STLEZ, you will be using the values - summarized below. See documentation for STL for definitions of - the parameters. - nt = smallest odd integer greater than or equal to - (1.5*np) / (1-(1.5/ns)). - nl = smallest odd integer greater than or equal to np. - nsjump = ns/10; ntjump = nt/10; nljump = nl/10. - ni = 2 if robust=.false.; else ni=1. - no = 0 if robust=.false.; else robustness iterations are carried - out until convergence of both seasonal and trend components, with - 15 iterations maximum. Convergence occurs if the maximum changes - in individual seasonal and trend fits are less than 1% of the - component's range after the previous iteration. //GO.SYSIN DD stlez.doc