Routine Methode.Autoregressive giving N, X
'' Calcule un intervalle a 0.95 end
'' utilisant la methode de regression

   Define i, j, k, N, No, p                as integer variable
   Define Critical.Value, Demie.Largeur, T as real    variable
   Define Df, V, W, Xbar                   as real    variable
   Define bsum, Lz, R, S, Vxbar, X         as a 1-dim real array
   Define b                                as a real 2-dim array


   if N < 50,
      Print 1 line thus
    Erreur ---- nombre d'observations insuffisant
      return
   endif

   Reserve S (*)                          as 51
   Reserve b(*,*)                         as 51 by *,
           bsum(*), Lz(*), R(*), Vxbar(*) as 51


   for i = 1 to N,
   do
      compute Xbar as the mean     of X(i)
      compute R(1) as the variance of X(i)
   loop

   let S(1) = R(1)

   skip 5 lines
   print 5 lines with TIME.V, Xbar thus
                 Temps de simulation *******.*****

          Estimation d'un intervalle a 0.95 --- Approche autoregressive

               Moyenne echantillonnee = ********.*****
   skip 3 lines
   begin report

    begin heading

         Print 3 lines thus
          Variance echan.    Intervalle (0.95)    Degree            Valeur
  Ordre   de la moy. echan.  Inf.       Sup.      liberte   T       Critique
  ------  -----------------  -----------------    --------  --     ---------
    end ''heading

      Reserve b(1, *) as 1
      let b(1, 1) = 1.0

      for i = 1 to 50,
      do
         let k = i + 1
         reserve b( k, *) as k

         for j = 1 to N-i,
            compute  R(k) as the mean of (X(j) - Xbar) *(X(i+j) - Xbar)

         for j = 1 to i,
         do
            compute V as the sum of b(i, j) * R(j)
            compute W as the sum of b(i, j) * R(i-j+2)
         loop

         let b(k, k) = - W / V
         let b(k, 1) = 1.0

         if i > 1,
            for j = 2 to i,
               let b(k, j) = b(i, j) + b(k, k) * b(i, k-j+1)
         endif

         for j = 1 to k,
         do
            compute S(k) as the sum of  b(k, j) * R(j)
            if j < k,
               compute bsum(i) as the sum of b(i, j)
               compute Lz(i)   as the sum of
                              (i - 1.0 - 2.0 * (j - 1.0)) * b(i, j)
            endif
         loop
         let Vxbar(i) = S(i) / (N * bsum(i) * bsum(i))
      loop

      let Df = N - 1.0

      for j = 1 to 50
      do'' La statistique T converge vers une Chi-carree
         print 1 line with s(51), j, s(j) thus
   s(51) = *******.***********   s(***) = *********.************
         Let T = N * (1.0 - S(51) / S(j))
         let Critical.Value =  ChiCarree.f ( REAL.F ( 51 - j), 0.90)
         if j > 1,
            Let Df = MAX.F ( 1.0, N * bsum(j) / (2.0 * Lz(j)))
         endif
         Let Demie.Largeur = Student.F ( Df) * Sqrt.F ( Vxbar(j))
         print 1 line with j-1, Vxbar(j),
                           Xbar-Demie.Largeur,  Xbar + Demie.Largeur,
                           Df, T, Critical.Value Thus
 ****  ********.******** *****.**** *****.****    **** ****.**   ****.**
         skip 1 line
         if T <= Critical.Value AND No = 0,
            let p  = j
            Let No = 1
         endif
      loop

      if p = 0,
         Print 1 line thus
          Aucun ordre ne satisfait le test
      else
         Print 1 line  with p-1 thus
          Ordre = *** satisfait le test
      endif
   end '' report

   skip 5 lines
   print 1 line thus
               Autocorrelation echantillonnee
   skip 1 line
   for i = 1 to p, write R(i) / R(1) as (4) E(18,6)
   skip 3 lines
   print 1 line thus
               Echantillon coefficient de correlation
   for i = 1 to p, write b(p, i) as (4) E(18,6)
   if p = 1,
      Release b(*,*), bsum(*), Lz(*), R(*), S(*), Vxbar(*)
      return
   endif


   Define Y as a 1-dim real array
   Define Theta, Z, ZZ as real variable

   Reserve Y(*) as p-1
   for i = 1 to p-1, let Y(i) = R(i) / R(1)

   let j = 1
   let Theta = 0.05
   let ZZ = RINF.C
   while (ZZ > 0.01)
   do
      let Z = 0
      for i = 1 to p-1,
      do
         subtract b(p, p+1-i) * Y(i) from Z
         if i < p-1,
            let y(i) = Y(i+1)
         endif
      loop

      let Y (p-1) =  Z
      let ZZ = ABS.F (Z)
      if ZZ <= Theta,
         skip 2 lines
         print 1 line with j, theta thus
         J minimum = ***** pour ABS.F(autocorrelation(j)) <= *.**
         if ZZ > 0.01,
            let Theta = 0.01
         endif
      endif
      Add 1 to j
   loop

   release bsum(*), Lz(*),R(*), S(*), Vxbar(*)

   for i = 1 to p-1,
      let Y(i) = -1.0
   let j = 1
   let Theta = 0.05
   let ZZ = RINF.C

   while (ZZ > 0.01)
   do
      let Z = 0
      for i = 1 to p-1
      do
         subtract b(p, p+1-i) * Y(i) from Z
         if i < p-1,
            let y(i) = Y(i+1)
         endif
      loop

      let Y (p-1) =  Z
      let ZZ = ABS.F (Z)
      if ZZ <= Theta,
         skip 2 lines
         print 1 line with j, theta thus
         J minimum = ***** pour ABS mean(j) cond - mean) / mean <= *.***
         if ZZ > 0.01,
            let Theta = 0.01
         endif
      endif
      Add 1 to j
   loop
   release b(*,*), Y(*)
   return
end '' Methode.Autoregressive

FUNCTION ChiCarree.f ( Df, P )

''*****************************************************************************
''* Fonction: Fonction qui calcule le P i
me quantile d'une distribution      *
''*           Chi carre avec Df degr de libert                             *
''* Usage   :                                                                 *
''* Entre  : Df : nombre de degr de libert                                 *
''*           P  : quantile   0.0 <= P <= 1.0                                 *
''* Retour  : Quantile                                                        *
''* E/S     : Aucune                                                          *
''* Portable: SIMSCRIPT II.5  (IBM 370 & IBM PC (Dos) & SUN)                  *
''* Notes   : Extrait de Principles of Discrete Event Simulation p. 254       *
''* Statut  : Dveloppement                                                   *
''* Creation:   05-92   Granger Louis                                         *
''*****************************************************************************

   DEFINE i                                       AS AN INTEGER VARIABLE
   DEFINE A, Df, P, Q, T, Xp, Num, Den,
          Y, Sqdf, sqhalf, Y2, Y3, Y4, Y5, Y6, Y7 AS REAL VARIABLES
   DEFINE C, H, D                                 AS 1-DIM REAL ARRAYS


   RESERVE C(*), D(*) AS 3,
           H(*)       AS 7

   LET Num = 0
   LET Den = 1.0
   LET Q   = P
   IF P > 0.5,
      LET Q = 1.0 - P
   ENDIF

   LET T = SQRT.F ( LOG.E.F ( 1.0 / (Q*Q)))

   LET C(1) = 2.515517
   LET C(2) = 0.802853
   LET C(3) = 0.010328

   LET D(1) = 1.432788
   LET D(2) = 0.189269
   LET D(3) = 0.001308

   FOR i = 1 TO 3
   DO
      ADD C(i) * T ** ( i-1 ) TO Num
      ADD D(i) * T ** i TO Den
   LOOP

   LET Xp = T - Num / Den
   IF P < 0.5,
      LET Xp = - Xp
   ENDIF

   LET Y      = Xp
   LET Sqdf   = SQRT.F ( Df )
   LET Sqhalf = SQRT.F ( 0.5 )
   LET Y2     = Y * Y
   LET Y3     = Y2 * Y
   LET Y4     = Y3 * Y
   LET Y5     = Y4 * Y
   LET Y6     = Y5 * Y
   LET Y7     = Y6 * Y


   LET H(1) = Y / Sqhalf
   LET H(2) = 2.0 * ( Y2 - 1.0 ) / 3.0
   LET H(3) = ( Y3 - 7.0 * Y ) * Sqhalf / 9.0
   LET H(4) = - ( 6.0 * Y4 + 14.0 * Y2 - 32.0 ) / 405.0
   LET H(5) = ( 9.0 * Y5 + 256.0 * Y3 - 433.0 * Y )
              * Sqhalf / 4860.0
   LET H(6) = ( 12.0 * Y6 - 243.0 * Y4 - 923.0 * Y2 + 1472.0 )
              / 25515.0
   LET H(7) =  - ( 3753.0 * Y7 + 4353.0 * Y5  - 289517.0 * Y3 -
               289717.0 * Y)  * Sqhalf / 9185400.0
   FOR i = 1 TO 7 COMPUTE A AS THE SUM OF H(i) / Sqdf ** i

   RETURN WITH ( A + 1.0) * Df

END '' ChiCarree





Routine Methode.Spectrale giving N, X
'' Cette routine utilise methode spectrale (Tukey-Hanning)
'' pour calculer un intervalle de confiance a 0.95

   Define  i, j, k, N                     as integer variable
   Define Df, Demie.Largeur, V, Vxbar, Xbar as real variable
   Define R, X as 1-dim real arrays

   if N < 512,
      print 1 line thus
      Erreur -- Nombre insuffisant d'observations
      return
   endif

   Reserve R(*) as 513

   for i = 1 to N,
      compute Xbar as the mean of X(i)

   skip 2 lines
   print 5 lines with TIME.V, Xbar thus
               Temps de simulation   *******.******

          Estimation d'un intervalle a 0.95 --- Approche spectrale

               Moyenne echantillonnee = ********.*****
   skip 2 lines
   begin report
    begin heading
         Print 3 lines thus
  No.     Degre    Variance echan.    Intervalle (0.95)
  lags    liberte  de la moy. echan.  Inf.       Sup.
  ------  -------  -----------------  -----------------
    end  ''heading
      for j = 1 to 513,
      do
         for i = 1 to N-j+1,
            compute R(j) as the sum of (X(i) - Xbar)*(X(i+j-1)-Xbar)
         let R(j) = R(j) / N
      loop
      let k = 1
      for i = 1 to 9,
      do
         let k = 2 * k
         for j = 1 to k,
            compute V as the sum of ( 1.0 + COS.F(J * PI.C / k)) * R(j+1)
         let Vxbar = ( V + R(1)) / ( N - k )
         let Df = 2.667 * N / k
         let Demie.Largeur = Student.F (Df) * SQRT.F  (Vxbar)
         skip 1 line
         print 1 line with  k, Df, Vxbar,
                            Xbar - Demie.Largeur,  Xbar + Demie.Largeur
                       thus
  ******* *******  *******.*********** ***.****** ***.******
      loop
      release R
      return
   end '' report
end '' Methode.Spectrale