# to unbundle, sh this file (in an empty directory)
echo Makefile 1>&2
-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
-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
-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
-.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
- 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
- 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
- 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
-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<nsh) {
- nleft = 1
- nright = len
- }
- else 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(maxs<work(i,6)) maxs = work(i,6)
- if(maxt<work(i,7)) maxt = work(i,7)
- if(mins>work(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(maxds<difs) maxds = difs
- if(maxdt<dift) maxdt = dift
- }
- if((maxds/(maxs-mins)<.01) & (maxdt/(maxt-mint)<.01)) break
- }
- }
-if(!robust){
- do i = 1,n
- rw(i) = 1.0
- }
-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 (n<0||ni<0) return #ERROR(number of items<0)
-if(n<2 | ni==0) return
-jl = 1
-ju = ni
-indl(1) = 1
-indu(1) = ni
-# arrays indl, indu keep account of the portion of ind related to the
-# current segment of data being ordered.
-i = 1
-j = n
-m = 1
-repeat {
- if (i<j) go to 10
- repeat {
- m = m-1
- if (m==0) break 2
- i = il(m)
- j = iu(m)
- jl = indl(m)
- ju = indu(m)
- if (jl<=ju) {
- while (j-i>10) {
-# 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(j)
- a(j) = t
- t = a(ij)
- if (a(i)>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
-.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
- 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
.
