[CONTACT]

[ABOUT]

[POLICY]

echo Makefile out out SYSIN DD

Found at: ftp.icm.edu.pl:70/packages/netlib/a/stl

# 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

.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]