README*README*README*README*README*README*README*README*README*README

There are 4 files here which should be separated:

	(1) gnpgnp.dat
	(2) D.par
	(3) gnpgnpD.out
	(4) arfima.f

(1) US GNP data used in Sowell (1992) and Koop et al. (1995).
(2) An example [AR(1)] of par file for arfima.exe
(3) The [arfima.exe + D.par] output file
(4) The FORTRAN code used in Koop et al. (1995).
    Note:  See also Chung (1994) for more efficient numerical
    methods.
    
Refs:    
    
C.-F.Chung (1994): "A Note on Calculating the Autocovariances
	of the Fractionally Intergrated ARMA Models,"
	Economics Letters 45, 293-297.   

G.Koop, E.Ley, J.Osiewalski and M.F.J.Steel (1995):
 	"Bayesian Analysis of Long Memory and Persistence using
	ARFIMA Models," Journal of Econometrics, forthcoming.
	http://econwpa.wustl.edu/listings/em//9505001

F.Sowell (1992): "Maximum Likelihood Estimation  of Stationary 
	Univariate Fractionally Integrated Time Series Models", 
	Journal of Econometrics, 53 (1992), 165-188. 
	
Email:
	gkoop@epas.utoronto.ca
	edley@eco.uc3m.es
	eeosiewa@cyf-kr.edu.pl
	steel@kub.nl

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ---> gnpgnp.dat
0.63216703E-02
0.36614602E-02
0.12016352E-01
0.62765574E-02
0.17602716E-01
0.91793301E-02
0.82077363E-02
-.11707538E-01
-.58603603E-02
0.75669220E-02
-.99207163E-02
0.39889810E-01
0.28166038E-01
0.36824175E-01
0.28084727E-01
0.20732641E-01
0.25931641E-01
0.22021712E-01
0.45827562E-02
0.96874382E-02
-.24130756E-02
0.89647479E-02
0.20539322E-01
0.17335365E-01
0.93871287E-02
-.46477833E-02
-.80983408E-02
-.13976369E-01
-.39886093E-02
0.11918416E-01
0.14562005E-01
0.21180822E-01
0.10895787E-01
0.13239027E-01
0.87296368E-02
-.19773273E-02
0.45420214E-02
0.72218760E-03
0.11030363E-01
0.82097491E-02
-.57957950E-03
0.58447772E-02
-.15619267E-01
-.20504103E-01
0.53637185E-02
0.23367684E-01
0.23401457E-01
0.12339263E-01
0.18869648E-01
-.45920791E-02
0.84940472E-02
0.17013985E-01
-.28756310E-02
0.95946277E-03
-.86080291E-02
0.10344712E-01
0.12368594E-01
0.14200450E-01
0.22241064E-01
0.13021017E-01
0.10351770E-01
0.92534252E-02
-.16559951E-02
0.13444382E-01
0.13750014E-01
0.17322219E-01
0.71605634E-02
0.22103214E-01
0.85333034E-02
0.10023878E-01
0.42725441E-02
0.21436835E-01
0.14378919E-01
0.15795993E-01
0.22746984E-01
0.19596265E-01
0.25992399E-02
0.10194692E-01
0.49016399E-02
0.56363379E-02
0.59595462E-02
0.14308285E-01
0.56230124E-02
0.11538839E-01
0.16872285E-01
0.77438205E-02
-.96470448E-03
0.13960014E-01
0.13646798E-02
0.55223713E-02
-.39944872E-02
-.61671105E-02
-.87225608E-03
0.12101835E-01
-.90729755E-02
0.26491615E-01
-.80693968E-04
0.51111893E-02
-.40143715E-04
0.21682142E-01
0.19258673E-01
0.10350472E-01
0.18589722E-01
0.23200493E-01
0.25570790E-02
-.98552743E-03
0.89073683E-02
-.55896483E-02
0.28350253E-02
-.13115541E-01
-.88278779E-02
-.19745495E-01
0.10127527E-01
0.16826472E-01
0.13827128E-01
0.18607364E-01
0.44473772E-02
0.41449001E-02
0.99202277E-02
0.13628357E-01
0.15997053E-01
0.19884581E-01
-.25684232E-02
0.87786950E-02
0.31095654E-01
0.85324477E-02
0.12333732E-01
0.31430233E-04
-.94333697E-03
0.89883319E-02
-.19036628E-02
0.99772379E-02
-.23912006E-01
0.66496733E-03
0.12613602E-01
0.19163783E-01
-.33480289E-02
0.44207108E-02
-.14066491E-01
-.15212989E-01
0.29919869E-02
-.80197449E-02
0.15204792E-02
0.86040336E-02
0.22251069E-01
0.14654397E-01
0.17597696E-01
0.25409183E-01
0.13324513E-01
0.64400530E-02
0.41384452E-02
0.11894358E-01
0.60473661E-02
0.10117571E-01
0.72895009E-02
0.15900653E-01
-.44440328E-02
0.21032770E-02
0.56943479E-02
0.13144430E-01
0.10648888E-01
0.12811526E-01
0.16085588E-01
0.99110844E-02
0.89913573E-02
0.79469959E-02
0.65827829E-02
0.91485680E-02
0.62384143E-02
0.73293961E-02
0.22794630E-02
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ---> D.par
D
-917364                          IDUM
25000                          NUMREP
330.0d0   
gnpgnp
171                            MYNOBS
1  0                           P AND Q
.FALSE.                        TRUE IF D FIXED
0.0D0              D       
-1.0d0  0.5d0   60             D
-1.0d0  0.5d0   30              PHI1
0.0d0 4.0d0 60              I(4)
0.0d0 6.5d0 60              I(12)
0.0d0 7.5d0 60              I(40)
2000.0d0 11000.0d0 60       T+40
16.0d0 20.0d0 60        SSE
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ---> gnpgnpD.out
 **********************************************************************
 **********************************************************************
          Model: D    ARFIMA(1,d,0)
          Data taken from: gnpgnp.dat
 **********************************************************************
 **********************************************************************
 ......................................................................
   random seed.............................. -917364
   number of antithetic replications........   25000
 ......................................................................
   order of ar..............................       1
   order of ma..............................       0
 ......................................................................
   total obs on differenced series..........     171
   mean of differenced series...............  0.0080
   standard deviation of diff series........  0.0107
 ......................................................................
   ff1 (parameter) .........................    330.0000
 ----------------------------------------------------------------------
   Number of times FRAC returns ierr.ne.0...      42
 ----------------------------------------------------------------------
   Accepted drawings in routine draw........100.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%               SSE*1000                         %%
              %%                                                %%
              %%   Mean...    16.80        STD...0.2837         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                              XXXX|   16.500 0.0185   1.8496
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   16.583 0.2199  23.8441
       XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   16.667 0.1895  42.7961
                  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   16.750 0.1421  57.0013
                            XXXXXXXXXXXXXXXXXXXXXX|   16.833 0.0986  66.8575
                                 XXXXXXXXXXXXXXXXX|   16.917 0.0759  74.4450
                                     XXXXXXXXXXXXX|   17.000 0.0584  80.2883
                                       XXXXXXXXXXX|   17.083 0.0498  85.2649
                                          XXXXXXXX|   17.167 0.0373  88.9996
                                           XXXXXXX|   17.250 0.0300  92.0033
                                             XXXXX|   17.333 0.0212  94.1206
                                              XXXX|   17.417 0.0166  95.7807
                                               XXX|   17.500 0.0121  96.9937
                                                XX|   17.583 0.0084  97.8311
                                                 X|   17.667 0.0062  98.4506
                                                 X|   17.750 0.0053  98.9822
                                                 X|   17.833 0.0032  99.2988
                                                 X|   17.917 0.0023  99.5266
                                                  |   18.000 0.0017  99.6919
                                                  |   18.083 0.0012  99.8100
                                                  |   18.167 0.0007  99.8822
                                                  |   18.250 0.0005  99.9272

              Mass to the left of           15.0000 is.....  0.0000 %
              Mass to the right of          20.0000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%             GNP_{T+40}                         %%
              %%                                                %%
              %%   Mean...    5790.        STD... 558.8         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                                  | 3800.000 0.0004   0.1007
                                                  | 3950.000 0.0006   0.1638
                                                  | 4100.000 0.0010   0.2679
                                                 X| 4250.000 0.0017   0.4421
                                                 X| 4400.000 0.0030   0.7386
                                                XX| 4550.000 0.0051   1.2507
                                              XXXX| 4700.000 0.0089   2.1408
                                            XXXXXX| 4850.000 0.0154   3.6777
                                       XXXXXXXXXXX| 5000.000 0.0259   6.2664
                                 XXXXXXXXXXXXXXXXX| 5150.000 0.0417  10.4341
                        XXXXXXXXXXXXXXXXXXXXXXXXXX| 5300.000 0.0628  16.7132
              XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 5450.000 0.0868  25.3958
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 5600.000 0.1085  36.2425
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 5750.000 0.1212  48.3610
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 5900.000 0.1208  60.4368
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 6050.000 0.1079  71.2238
              XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX| 6200.000 0.0874  79.9649
                       XXXXXXXXXXXXXXXXXXXXXXXXXXX| 6350.000 0.0653  86.4974
                               XXXXXXXXXXXXXXXXXXX| 6500.000 0.0459  91.0836
                                     XXXXXXXXXXXXX| 6650.000 0.0308  94.1638
                                          XXXXXXXX| 6800.000 0.0201  96.1759
                                             XXXXX| 6950.000 0.0130  97.4720
                                               XXX| 7100.000 0.0083  98.3048
                                                XX| 7250.000 0.0054  98.8430
                                                 X| 7400.000 0.0035  99.1951
                                                 X| 7550.000 0.0023  99.4292
                                                 X| 7700.000 0.0016  99.5877
                                                  | 7850.000 0.0011  99.6971
                                                  | 8000.000 0.0008  99.7740
                                                  | 8150.000 0.0006  99.8292
                                                  | 8300.000 0.0004  99.8693
                                                  | 8450.000 0.0003  99.8991
                                                  | 8600.000 0.0002  99.9214

              Mass to the left of         2000.0000 is.....  0.0000 %
              Mass to the right of       11000.0000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%                delta                           %%
              %%                                                %%
              %%   Mean...  -0.2897        STD...0.2166         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                                 X|   -0.750 0.0010   0.1800
                                                XX|   -0.725 0.0020   0.3795
                                               XXX|   -0.700 0.0038   0.7639
                                              XXXX|   -0.675 0.0047   1.2376
                                           XXXXXXX|   -0.650 0.0084   2.0779
                                          XXXXXXXX|   -0.625 0.0089   2.9651
                                  XXXXXXXXXXXXXXXX|   -0.600 0.0176   4.7238
                            XXXXXXXXXXXXXXXXXXXXXX|   -0.575 0.0243   7.1509
                               XXXXXXXXXXXXXXXXXXX|   -0.550 0.0211   9.2573
                       XXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.525 0.0299  12.2482
                XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.500 0.0377  16.0135
           XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.475 0.0435  20.3640
                XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.450 0.0386  24.2209
           XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.425 0.0438  28.5964
         XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.400 0.0455  33.1505
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.375 0.0503  38.1833
  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.350 0.0544  43.6235
           XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.325 0.0436  47.9877
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.300 0.0562  53.6053
         XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.275 0.0462  58.2254
           XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.250 0.0438  62.6074
      XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.225 0.0497  67.5787
                XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.200 0.0387  71.4483
                            XXXXXXXXXXXXXXXXXXXXXX|   -0.175 0.0249  73.9380
                        XXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.150 0.0293  76.8684
                            XXXXXXXXXXXXXXXXXXXXXX|   -0.125 0.0252  79.3924
                                XXXXXXXXXXXXXXXXXX|   -0.100 0.0203  81.4184
                               XXXXXXXXXXXXXXXXXXX|   -0.075 0.0208  83.4987
                                  XXXXXXXXXXXXXXXX|   -0.050 0.0182  85.3230
                                  XXXXXXXXXXXXXXXX|   -0.025 0.0181  87.1338
                                     XXXXXXXXXXXXX|    0.000 0.0150  88.6314
                                    XXXXXXXXXXXXXX|    0.025 0.0162  90.2501
                                     XXXXXXXXXXXXX|    0.050 0.0151  91.7568
                                     XXXXXXXXXXXXX|    0.075 0.0142  93.1740
                                        XXXXXXXXXX|    0.100 0.0111  94.2868
                                         XXXXXXXXX|    0.125 0.0102  95.3110
                                           XXXXXXX|    0.150 0.0075  96.0567
                                           XXXXXXX|    0.175 0.0073  96.7910
                                             XXXXX|    0.200 0.0058  97.3746
                                             XXXXX|    0.225 0.0051  97.8887
                                              XXXX|    0.250 0.0048  98.3713
                                               XXX|    0.275 0.0031  98.6812
                                               XXX|    0.300 0.0029  98.9724
                                                XX|    0.325 0.0027  99.2402
                                                XX|    0.350 0.0022  99.4624
                                                XX|    0.375 0.0017  99.6357
                                                 X|    0.400 0.0010  99.7360
                                                 X|    0.425 0.0009  99.8259
                                                 X|    0.450 0.0008  99.9029

              Mass to the left of           -1.0000 is.....  0.0000 %
              Mass to the right of           0.5000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%                Phi1                            %%
              %%                                                %%
              %%   Mean...  -0.6390        STD...0.1986         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                               XXX|   -0.950 0.0064   0.6417
                                    XXXXXXXXXXXXXX|   -0.900 0.0323   3.8682
                    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.850 0.0690  10.7717
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.800 0.1155  22.3189
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.750 0.1159  33.9080
   XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.700 0.1096  44.8642
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.650 0.1158  56.4437
       XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.600 0.0998  66.4254
                     XXXXXXXXXXXXXXXXXXXXXXXXXXXXX|   -0.550 0.0674  73.1615
                           XXXXXXXXXXXXXXXXXXXXXXX|   -0.500 0.0537  78.5339
                              XXXXXXXXXXXXXXXXXXXX|   -0.450 0.0474  83.2754
                                   XXXXXXXXXXXXXXX|   -0.400 0.0355  86.8209
                                      XXXXXXXXXXXX|   -0.350 0.0273  89.5485
                                     XXXXXXXXXXXXX|   -0.300 0.0311  92.6591
                                         XXXXXXXXX|   -0.250 0.0220  94.8610
                                           XXXXXXX|   -0.200 0.0156  96.4229
                                             XXXXX|   -0.150 0.0117  97.5972
                                              XXXX|   -0.100 0.0084  98.4362
                                               XXX|   -0.050 0.0065  99.0858
                                                XX|    0.000 0.0041  99.4912
                                                 X|    0.050 0.0027  99.7641
                                                 X|    0.100 0.0014  99.9044

              Mass to the left of           -1.0000 is.....  0.0000 %
              Mass to the right of           0.5000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%                I(4)                            %%
              %%                                                %%
              %%   Mean...    1.548        STD...0.2107         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                                 X|    1.067 0.0024   0.2911
                                               XXX|    1.133 0.0080   1.0953
                                          XXXXXXXX|    1.200 0.0218   3.2724
                                   XXXXXXXXXXXXXXX|    1.267 0.0397   7.2415
                    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.333 0.0784  15.0784
            XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.400 0.0993  25.0066
  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.467 0.1262  37.6261
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.533 0.1309  50.7166
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.600 0.1188  62.5966
        XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.667 0.1106  73.6566
                 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.733 0.0853  82.1819
                           XXXXXXXXXXXXXXXXXXXXXXX|    1.800 0.0607  88.2547
                                  XXXXXXXXXXXXXXXX|    1.867 0.0428  92.5300
                                       XXXXXXXXXXX|    1.933 0.0289  95.4156
                                           XXXXXXX|    2.000 0.0181  97.2216
                                              XXXX|    2.067 0.0113  98.3516
                                               XXX|    2.133 0.0072  99.0670
                                                XX|    2.200 0.0042  99.4855
                                                 X|    2.267 0.0023  99.7169
                                                 X|    2.333 0.0014  99.8527
                                                  |    2.400 0.0007  99.9248

              Mass to the left of            0.0000 is.....  0.0000 %
              Mass to the right of           4.0000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%                I(12)                           %%
              %%                                                %%
              %%   Mean...    1.323        STD...0.4055         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                                 X|    0.650 0.0041   0.4467
                                           XXXXXXX|    0.758 0.0193   2.3770
                                XXXXXXXXXXXXXXXXXX|    0.867 0.0498   7.3566
                 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    0.975 0.0919  16.5508
    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.083 0.1289  29.4389
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.192 0.1390  43.3342
    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.300 0.1292  56.2494
         XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.408 0.1150  67.7468
                    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.517 0.0840  76.1483
                             XXXXXXXXXXXXXXXXXXXXX|    1.625 0.0576  81.9062
                                 XXXXXXXXXXXXXXXXX|    1.733 0.0472  86.6263
                                     XXXXXXXXXXXXX|    1.842 0.0357  90.1925
                                         XXXXXXXXX|    1.950 0.0255  92.7422
                                            XXXXXX|    2.058 0.0176  94.5004
                                             XXXXX|    2.167 0.0134  95.8401
                                               XXX|    2.275 0.0094  96.7756
                                               XXX|    2.383 0.0089  97.6678
                                                XX|    2.492 0.0054  98.2081
                                                XX|    2.600 0.0046  98.6729
                                                 X|    2.708 0.0031  98.9817
                                                 X|    2.817 0.0026  99.2440
                                                 X|    2.925 0.0021  99.4529
                                                  |    3.033 0.0014  99.5899
                                                  |    3.142 0.0011  99.7001
                                                  |    3.250 0.0008  99.7804
                                                  |    3.358 0.0006  99.8415
                                                  |    3.467 0.0005  99.8903
                                                  |    3.575 0.0004  99.9283

              Mass to the left of            0.0000 is.....  0.0000 %
              Mass to the right of           6.5000 is.....  0.0000 %
 ----------------------------------------------------------------------

              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              %%                I(40)                           %%
              %%                                                %%
              %%   Mean...   0.9956        STD...0.6738         %%
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                                      UBound    PDF   CDF(%)
                                                  +-------------------------
                                                 X|    0.250 0.0039   0.3900
                                     XXXXXXXXXXXXX|    0.375 0.0378   4.1750
                 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    0.500 0.1002  14.1922
 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    0.625 0.1480  28.9914
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    0.750 0.1498  43.9745
      XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    0.875 0.1321  57.1854
                XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX|    1.000 0.1009  67.2798
                           XXXXXXXXXXXXXXXXXXXXXXX|    1.125 0.0703  74.3099
                                  XXXXXXXXXXXXXXXX|    1.250 0.0486  79.1728
                                      XXXXXXXXXXXX|    1.375 0.0374  82.9149
                                         XXXXXXXXX|    1.500 0.0255  85.4639
                                          XXXXXXXX|    1.625 0.0237  87.8341
                                           XXXXXXX|    1.750 0.0204  89.8775
                                             XXXXX|    1.875 0.0156  91.4361
                                              XXXX|    2.000 0.0127  92.7048
                                               XXX|    2.125 0.0104  93.7455
                                               XXX|    2.250 0.0102  94.7634
                                               XXX|    2.375 0.0076  95.5276
                                                XX|    2.500 0.0051  96.0336
                                                XX|    2.625 0.0055  96.5882
                                                XX|    2.750 0.0046  97.0514
                                                 X|    2.875 0.0034  97.3932
                                                 X|    3.000 0.0037  97.7598
                                                 X|    3.125 0.0030  98.0597
                                                 X|    3.250 0.0027  98.3332
                                                 X|    3.375 0.0019  98.5199
                                                 X|    3.500 0.0018  98.6999
                                                 X|    3.625 0.0017  98.8660
                                                  |    3.750 0.0012  98.9875
                                                  |    3.875 0.0011  99.0969
                                                  |    4.000 0.0011  99.2095
                                                  |    4.125 0.0011  99.3155
                                                  |    4.250 0.0010  99.4115
                                                  |    4.375 0.0009  99.4972
                                                  |    4.500 0.0006  99.5548
                                                  |    4.625 0.0005  99.6065
                                                  |    4.750 0.0006  99.6713
                                                  |    4.875 0.0004  99.7089
                                                  |    5.000 0.0003  99.7419
                                                  |    5.125 0.0002  99.7657
                                                  |    5.250 0.0004  99.8008
                                                  |    5.375 0.0003  99.8268
                                                  |    5.500 0.0002  99.8486
                                                  |    5.625 0.0002  99.8706
                                                  |    5.750 0.0002  99.8955
                                                  |    5.875 0.0002  99.9176

              Mass to the left of            0.0000 is.....  0.0000 %
              Mass to the right of           7.5000 is.....  0.0009 %
 ----------------------------------------------------------------------
>> sum of the weights:6460353420.6221
 ......................................................................
 >> moments of the weights:
   -----------------------------------------
            MEAN             STD
   -----------------------------------------
            0.1292E+06     0.5663E+06
 ----------------------------------------------------------------------
 >> posterior moments of beta:
   -----------------------------------------
         vname       MEAN      STD
   -----------------------------------------
       delta     -0.2897    0.2166
       Phi1      -0.6390    0.1986
 ----------------------------------------------------------------------
 >> posterior moments of impulse responses:
 -----------------------------------------------------------------
 I(i)         MEAN        STD      I(i)      MEAN          STD
 -----------------------------------------------------------------
 I( 1)      1.3492      0.0684  |  I(21)      1.1499      0.5250
 I( 2)      1.4916      0.1241  |  I(22)      1.1367      0.5353
 I( 3)      1.5421      0.1708  |  I(23)      1.1244      0.5453
 I( 4)      1.5483      0.2107  |  I(24)      1.1129      0.5549
 I( 5)      1.5324      0.2452  |  I(25)      1.1021      0.5642
 I( 6)      1.5060      0.2755  |  I(26)      1.0920      0.5731
 I( 7)      1.4749      0.3025  |  I(27)      1.0825      0.5817
 I( 8)      1.4424      0.3267  |  I(28)      1.0735      0.5901
 I( 9)      1.4103      0.3489  |  I(29)      1.0650      0.5981
 I(10)      1.3794      0.3692  |  I(30)      1.0569      0.6060
 I(11)      1.3502      0.3880  |  I(31)      1.0493      0.6136
 I(12)      1.3229      0.4055  |  I(32)      1.0421      0.6210
 I(13)      1.2974      0.4218  |  I(33)      1.0353      0.6282
 I(14)      1.2738      0.4372  |  I(34)      1.0288      0.6352
 I(15)      1.2519      0.4518  |  I(35)      1.0226      0.6420
 I(16)      1.2316      0.4655  |  I(36)      1.0166      0.6487
 I(17)      1.2127      0.4786  |  I(37)      1.0110      0.6552
 I(18)      1.1953      0.4910  |  I(38)      1.0056      0.6615
 I(19)      1.1790      0.5028  |  I(39)      1.0005      0.6677
 I(20)      1.1639      0.5141  |  I(40)      0.9956      0.6738
 ----------------------------------------------------------------------
 >> First 2 predictive moments of log(GNP_{T+i}):
   ------------------------------------------
   time               MEAN            STD
   ------------------------------------------
   T+ 1              8.3423         0.0100
   T+ 2              8.3496         0.0169
   T+ 3              8.3575         0.0226
   T+ 4              8.3656         0.0276
   T+ 5              8.3739         0.0319
   T+ 6              8.3822         0.0357
   T+ 7              8.3905         0.0391
   T+ 8              8.3989         0.0421
   T+ 9              8.4072         0.0449
   T+10              8.4155         0.0474
   T+11              8.4238         0.0498
   T+12              8.4321         0.0520
   ------------------------------------------
    (log(GNP_T) is:     8.3362 )
 ----------------------------------------------------------------------
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ---> arfima.f
c
c-----------------------------------------------------------------------
c ARFIMA.F: last modified: edley -> 24/April/95
c-----------------------------------------------------------------------
c
c  WARNING: All the variables must be initialized to 0.  Make sure
c           that you compile the code with an option that forces that.
c           (e.g., -K in the HP-AUX compiler.)
c
c This is ARFIMA.F by G.Koop, E.Ley, J.Osiewalski and M.Steel (1995):
c "Bayesian Analysis of Long Memory and Persistence using
c  ARFIMA Models," Journal of Econometrics, forthcoming.
c
c  gkopp@epas.utoronto.ca
c  edley@eco.uc3m.es
c  eeosiewa@cyf-kr.edu.pl
c  steel@kub.nl
c
c It is based on F.Sowell's code (see "Maximum Likelihood Estimation  
c of Stationary Univariate Fractionally Integrated Time Series Models", 
c Journal of Econometrics, 53 (1992), 165--188).  We thank F.Sowell
c for sharing his code with us.
c
c At present, it handles models up to order p=3 and q=3.  By modifying
c the subroutine DRAW (which draws the parameters) and, for d=0, 
c the FRAC routine (which evaluates the likelihood), any other model 
c can be used.
c
c-----------------------------------------------------------------------
c
c  - INPUT FILES:
c          		fdata............xxxxxx.dat
c          		fpara............arfima.par
c  - OUTPUT FILES:
c			fout.............arfimaX.out
c			ferr.............arfimaX.err
c			fmma.............arfimaX.mma
c
c  The execution is driven by fpara.  This file must contain:
c 
c   1. A character which will be appended to 'arfima' in the .out and 
c      .err files
c   2. A seed for the random number generator
c   3. Number of antithetic repetitions
c   4. ff1
c   5. A char*6 to identify the file to read the data in from 
c   6. The number of obs to be read from fdata
c   7. The order of the AR and MA polynomials
c   8. .false. if d is fixed; .true. otherwise
c   9. min, max and number of bins of the posterior plots
c      (one line for each parameter, I(i), etc; see below)
c
c  The strategy is to have many par files: A.par, B.par...
c  (each of them with letter A, B, etc, in the 1st position, 1st line)
c  and rename 'X.par' to 'arfima.par' when, say, model X 
c  wants to be run.  At runtime, the parameters will always be read
c  from arfima.par.  arfimaX.out and arfimaX.err will then contain the 
c  output.
c  
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      program arfima
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      implicit double precision (a-h,o-z)
      double precision impuls
      integer mynobs,p, q, np, ndraws, kdfix, nfails
      integer fdata,fpara,fout,ferr,fmma
      logical dfixed
c      
      character*1 letra, delta
      character*9 chnum
      character*6 nombre
c
      common/data/x(500)
      common/covar/sig(500),v2(500),chl(124750)
      common/parms/mynobs,p,q
      common/files/fdata,fpara,fout,ferr,fmma
      common/hat/uhat,sigtmp,sse,dm1
      common/seed/idum
c
c-----------------------------------------------------------------------
c  parameters:
c 	nimp...horizon for calculating impulse responses
c 	ndim...maximum dimension of np
c-----------------------------------------------------------------------
c
      integer nimp, ndim, maxnpred, df
      parameter(nimp=40, ndim=7, maxnpred=40)
      dimension betam(ndim), betasd(ndim),beta1(ndim),beta2(ndim)
      dimension din1(nimp+1),din2(nimp+1),impuls(nimp),sdimp(nimp)
c
c-----------------------------------------------------------------------
c  posterior densities stuff
c-----------------------------------------------------------------------
c
      parameter(yT=4172.4d0)
      integer maxnb, nbins, ndim3, ncdf
      parameter(rpmin=0.0d0,rpmax=4.0d0,nbrp=40)
      double precision scale, cdf, pdf, dummy, dat1, dat2
      double precision below, above, ave, ssq
      parameter (maxnb=100, ndim3=ndim+3, ncdf = ndim3*maxnb)
      dimension cdf(ndim+3,maxnb), pdf(maxnb)
      dimension dummy(maxnb)
      dimension dat1(3), dat2(3)
      dimension nbins(ndim3), dmax(ndim3), dmin(ndim3)
      dimension below(ndim3), above(ndim3), ave(ndim3), ssq(ndim3)     
      character*7 name(ndim3)
c ratio of polynomials
      dimension cdfrp(maxnb)
c predictives      
      dimension ypm(12),ypm1(13),ypm2(13)
      dimension ypsd(12),ypp1(13),ypp2(13)
      dimension cdfy40(maxnb)
c sse
      dimension cdfsse(maxnb)
c
c-----------------------------------------------------------------------
c  initialize stuff
c-----------------------------------------------------------------------
c
      data w,sdw/2*0.0d0/
      data ndraws,nd/2*0/
      data nfails/0/
c for hist and stats of parameters            
      data ((cdf(i,j), i=1,ndim3),j=1,maxnb)/ncdf*0.0d0/
      data (impuls(i), i=1,nimp)/nimp*0.0d0/ 
      data (sdimp(i), i=1,nimp)/nimp*0.0d0/ 
      data (betam(i), i=1,ndim)/ndim*0.0d0/
      data (betasd(i), i=1,ndim)/ndim*0.0d0/
      data (below(i), i=1,ndim3)/ndim3*0.0d0/ 
      data (above(i), i=1,ndim3)/ndim3*0.0d0/       
      data (ave(i), i=1,ndim3)/ndim3*0.0d0/       
      data (ssq(i), i=1,ndim3)/ndim3*0.0d0/
c ratio of polynomials
      data (cdfrp(i), i=1,maxnb)/maxnb*0.0d0/
      data brp,arp,averp,ssqrp/4*0.0d0/
c predictives
      data (ypm(i), i=1,12)/12*0.0d0/
      data (ypsd(i), i=1,12)/12*0.0d0/      
      data (cdfy40(i), i=1,maxnb)/maxnb*0.0d0/
      data by40,ay40,avey40,ssqy40/4*0.0d0/
c sse
      data (cdfsse(i), i=1,maxnb)/maxnb*0.0d0/
      data bsse,asse,avesse,ssqsse/4*0.0d0/
c
c-----------------------------------------------------------------------
c .par file
c-----------------------------------------------------------------------
c
       open(fpara,file='arfima.par')
c
c-----------------------------------------------------------------------
c read parameters:
c
c    letra.......character to be appended to 'nombre' in .out, .err
c    idum........random seed (negative integer)
c    numrep......number of antithetic replications (integer)
c    ff1.........
c    nombre......char*6 to identify data-input & output files
c    mynobs......total obs on differenced series (integer) 
c    p...........order of ar (integer) 
c    q...........order of ma (integer)
c    dfixed......logical (logical--true if d fixed; false otherwise)
c    dvalue......value of d (real; only if dfixed=.true.)
c
c-----------------------------------------------------------------------
c
      read(fpara,9) letra
    9 format(a1) 
      read(fpara,*) idum
      read(fpara,*) numrep
      scale = 10.0d0/(numrep*1.0d0)
      read(fpara,*) ff1
      read(fpara,11) nombre
   11 format(a6)
      read(fpara,*) mynobs
      df=mynobs-1
      read(fpara,*) p, q
      np = p + q + 1
      read(fpara,*) dfixed
c
      if (dfixed) then 
	kdfix=1
      else
        kdfix=0
      endif
      read(fpara,*) dvalue
c
c-----------------------------------------------------------------------
c  read min, max and number of bins for histogram plots.
c
c   (1) first for each of the components in: 
c   	beta = (d, ar-phi-vector, ma-theta-vector)
c   (2) then for I(4), I(12) and I(40)
c
c              min     max    nbins
c
c     d         (if dfixed = .false.)
c     phi1
c     ...
c     phip
c     theta1
c     ...
c     thetaq
c     i(4)
c     i(12)
c     i(40)
c
c-----------------------------------------------------------------------
c
      do 15 i=1,np+3
   15 read(fpara,*) dmin(i),dmax(i),nbins(i)
c
c predictives' stuff
c
      read(fpara,*) ymin40, ymax40, nbinsy
c
c SSE
c
      read(fpara,*) ssemin, ssemax, nbinsse
c
c  close the .par file (to be able to run various processes in
c  paralell)
c
      close(fpara)
c
c-----------------------------------------------------------------------
c  output files
c-----------------------------------------------------------------------
c
      open(fout,file=nombre//letra//'.out')
      open(ferr,file=nombre//letra//'.err')
      open(fmma,file=nombre//letra//'.mma')
c
c-----------------------------------------------------------------------
c  names of variables in post plots
c-----------------------------------------------------------------------
c
      name(1) = 'delta'
c 
      chnum='123456789'
c       
      do 20 i=1,p
   20 name(i+1) = 'Phi'//chnum(i:i)
c       
      do 21 i=1,q
   21 name(i+1+p) = 'Theta'//chnum(i:i)
c
      name(np+1) = 'I(4)'
      name(np+2) = 'I(12)'
      name(np+3) = 'I(40)'
c       
c-----------------------------------------------------------------------
c  read  data and compute its mean and standard variation
c-----------------------------------------------------------------------
c
      open(fdata,file=nombre//'.dat')
c
      sum = 0.0d0
      sq = 0.0d0
      do 40 i =1,mynobs
        read(fdata,*) x(i)
        sum = sum + x(i)
   40 sq = sq + x(i)**2
c
c  close the .dat file.
c
      close(fdata)
c
      dmean = sum/mynobs
      scmom = sq/mynobs
      vardat = scmom - dmean**2
      sddat = sqrt(vardat)
c
c-----------------------------------------------------------------------
c  write misc input info out to arfima.out
c-----------------------------------------------------------------------
c
      call barra('*',70,fout)
      call barra('*',70,fout)
      if (dfixed) then
      delta='1'
      else
      delta='d'
      endif
      write(fout,50) letra,p,delta,q
   50 format(10x,'Model: ',a1,4x,'ARFIMA(',i1,',',a1,',',i1,')')   
      write(fout,51) nombre//'.dat'
   51 format(10x,'Data taken from: ',a10)   
      call barra('*',70,fout)
      call barra('*',70,fout)
c
      if (dfixed) then
        call barra('.',70,fout)
        write(fout,55) dvalue
   55   format(3x,'conditional on delta equal to............',f8.4)
      endif
c
      call barra('.',70,fout) 
c
      write(fout, 110) idum
  110 format(3x,'random seed..............................',i8)
      write(fout, 120) numrep
  120 format(3x,'number of antithetic replications........',i8)
c
      call barra('.',70,fout) 
c
      write(fout, 130) p
  130 format(3x,'order of ar..............................',i8)
      write(fout, 140) q
  140 format(3x,'order of ma..............................',i8)
c
      call barra('.',70,fout) 
c
      write(fout, 200) mynobs
  200 format(3x,'total obs on differenced series..........',i8)
      write(fout, 210) dmean
  210 format(3x,'mean of differenced series...............',f8.4)
      write(fout, 220) sddat
  220 format(3x,'standard deviation of diff series........',f8.4) 
c
      call barra('.',70,fout) 
c
      write(fout, 230) ff1
  230 format(3x,'ff1 (parameter) .........................',f12.4)
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c ---------------   start the big monte carlo loop ---------------------
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c      seg = 0.0
c      seg = secnds(seg)
c
c Starting the Big MonteCarlo Loop'
c
      do 2000 iloop=1, numrep
c
c take two antithetic draws from beta call draws beta1 and beta2
c
  300 nd = ndraws
      call draw(beta1, beta2, dfixed, dvalue, p, q, nd, zlit1, zlit2)
c
c evaluate log of posterior at the draws
c
      call frac(beta1,f1,maxnpred,ierr)
      if (ierr.ne.0) then
          nfails = nfails+1
          goto 300
	else
	 do 310 i=1,12
 310     call predict(i,ypm1(i),ypp1(i))
         call predict(40,ypm1(13), ypp1(13))
	 sse1 = sse*1000
      endif
      call frac(beta2,f2,maxnpred,ierr)
      if (ierr.ne.0) then
          nfails = nfails+1
          goto 300
        else
 	 do 320 i=1,12
 320     call predict(i,ypm2(i),ypp2(i))
         call predict(40,ypm2(13), ypp2(13))
	 sse2 = sse*1000
      endif
c      
      ndraws = nd
c
c get impulse-response fcns
c
      call zimp(beta1,nimp,din1)
      call zimp(beta2,nimp,din2)
c
c calculate weight functions for monte carlo
c
      w1 = f1 - ff1 - zlit1
      w2 = f2 - ff1 - zlit2
c      
      if (w1.gt.20)  then
       write(ferr,*) '*** log(w1) bigger than 20!!!...', w1
       do 1000 i=1+kdfix,np
 1000  write(fout,*) name(i), beta1(i)
       write(fout,*)
      endif
c
      w1=exp(w1)
c      
      if (w2.gt.20) then
       write(ferr,*) '*** log(w2) bigger than 20!!!...', w2
       do 1010 i=1+kdfix,np
 1010  write(fout,*) name(i), beta2(i)
       write(fout,*)
      endif
c
      w2=exp(w2)
c
c-----------------------------------------------------------------------
c
c calculate means and standard deviations of the parameters
c & accumulate weights in array cdf for posterior plots
c
      do 1300 i=1+kdfix,np
c
c moments --- they're also computed by hist; 
c we leave these redundant statements for checking purposes   
c
      betam(i)=betam(i) + w1*beta1(i) + w2*beta2(i)
      betasd(i)=betasd(i) + w1*beta1(i)**2 + w2*beta2(i)**2
c
c cdf's      
c      
      do 1100 j=1, nbins(i)
 1100 dummy(j) = cdf(i,j)
      call acc(beta1(i),w1*scale, dmin(i), dmax(i), nbins(i), dummy,
     & below(i), above(i), ave(i), ssq(i))
      call acc(beta2(i),w2*scale, dmin(i), dmax(i), nbins(i), dummy,
     & below(i), above(i), ave(i), ssq(i))
      do 1200 j=1, nbins(i)
 1200 cdf(i,j) = dummy(j)
 1300 continue
c
c-----------------------------------------------------------------------
c
c Ratio of MA/AR polynomials at 1
c
      if ((dfixed).and.(dvalue.eq.0)) then
      call acc(ratpol(beta1,p,q),w1*scale,rpmin,rpmax,nbrp,cdfrp,
     & brp,arp,averp,ssqrp) 
      call acc(ratpol(beta2,p,q),w2*scale,rpmin,rpmax,nbrp,cdfrp,
     & brp,arp,averp,ssqrp) 
      endif
c
c-----------------------------------------------------------------------
c
c same for impulse-response fcns I(4), I(12), I(40)
c
      dat1(1) = din1(4)
      dat1(2) = din1(12)
      dat1(3) = din1(40)
      dat2(1) = din2(4)
      dat2(2) = din2(12)
      dat2(3) = din2(40)
c      
      do 1600 i=np+1,np+3
      do 1400 j=1,nbins(i)
 1400 dummy(j) = cdf(i,j)
      call acc(dat1(i-np),w1*scale, dmin(i), dmax(i), nbins(i), dummy,
     & below(i), above(i), ave(i), ssq(i))
      call acc(dat2(i-np),w2*scale, dmin(i), dmax(i), nbins(i), dummy,
     & below(i), above(i), ave(i), ssq(i))
      do 1500 j=1,nbins(i)
 1500 cdf(i,j) = dummy(j)
 1600 continue         
c
c-----------------------------------------------------------------------
c
c SSE
c
      call acc(sse1,w1*scale, ssemin, ssemax,
     & nbinsse, cdfsse,bsse,asse,avesse,ssqsse)
      call acc(sse2,w2*scale, ssemin, ssemax,
     & nbinsse, cdfsse,bsse,asse,avesse,ssqsse)
c-----------------------------------------------------------------------
c
c  2 first predictive moments for log(GNP_{T+i}); i=1,2,...,12
c
       do 1650 i=1,12
       ypm(i) = ypm(i) + ypm1(i)*w1 + ypm2(i)*w2
 1650  ypsd(i) = ypsd(i) + w1/sqrt(ypp1(i)) + w2/sqrt(ypp2(i))
c
c  f[y_{T+40}|data]
c
      step = (ymax40-ymin40)/nbinsy
      y = ymin40+step/2.0d0
      do 1700 j=1,nbinsy
      ww = w1*scale*studt(log(y),df,ypm1(13)+log(yT),ypp1(13))/y
      call acc(y,ww,ymin40,ymax40,nbinsy,cdfy40,by40,ay40,avey40,ssqy40)
      ww = w2*scale*studt(log(y),df,ypm2(13)+log(yT),ypp2(13))/y
      call acc(y,ww,ymin40,ymax40,nbinsy,cdfy40,by40,ay40,avey40,ssqy40)
 1700 y = y+step
c
c
c-----------------------------------------------------------------------
c
c  compute mean and sd of weights
c     
      w = w + w1 + w2
      sdw = sdw+w1**2+w2**2
c
c-----------------------------------------------------------------------
c
c calculate mean and sd of impulse response function
c
      do 1750 i=1,nimp
        impuls(i) = impuls(i) + w1*din1(i) + w2*din2(i)
        sdimp(i) = sdimp(i) + w1*din1(i)**2 + w2*din2(i)**2
 1750 continue
c
 2000 continue
c
c      seg = secnds(seg)
c      write(fout,*) 'End of Big MonteCarlo Loop'
c      write(fout,*) '*---> ', seg, ' seconds'
c      write(fout,*) '      (', seg/60.0, ' minutes)'
c      write(fout,*) '      (', seg/3600.0, ' hours)'
c      
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c ----------------   end of biggg monte carlo looop  -------------------
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
c  ierr's in FRAC and acceptance rate in DRAW
c
      call barra('-',70,fout)
      write(fout,2009) nfails
 2009 format(3x,'Number of times FRAC returns ierr.ne.0...',i8)
      call barra('-',70,fout)
      write(fout,2010) ((numrep*1.0d0)/(ndraws*1.0d0))*100.0d0
 2010 format(3x,'Accepted drawings in routine draw........',f8.4,' %')
      write(fout,*) '(Corrected code!)'
c
c  posterior plots: SSE
c
      call barra('-',70,fout)
      call hist('SSE*1000',letra, cdfsse, pdf, ssemin, ssemax, 
     & nbinsse, bsse, asse, avesse, ssqsse)
c
c  posterior plots: predictives  40
c
      call barra('-',70,fout)
      call hist('GNP_{T+40}',letra, cdfy40, pdf, ymin40, ymax40, 
     & nbinsy, by40, ay40, avey40, ssqy40)
c
c  posterior plots: parameters in beta and I(4), I(12), I(40)
c
      do 2020 i=1+kdfix,np+3
      do 2015 j=1,nbins(i)
 2015 dummy(j) = cdf(i,j)
      call barra('-',70,fout)
 2020 call hist(name(i), letra, dummy, pdf, dmin(i), dmax(i), nbins(i),
     & below(i), above(i), ave(i), ssq(i))
c
c  posterior plots: ratio of polynomials
c
      if ((dfixed).and.(dvalue.eq.0)) then
      call barra('-',70,fout)
      call hist('I(inf)', letra, cdfrp, pdf, rpmin,rpmax,nbrp,
     & brp,arp,averp,ssqrp)
      endif
c 
c-----------------------------------------------------------------------
c
c compute means and sd's of beta and I()'s
c
      do 2060 i=1+kdfix,np
       betam(i)=betam(i)/w
       betasd(i)=sqrt((betasd(i)/w) - betam(i)**2)
 2060 continue
      do 2070 i=1,nimp
        impuls(i)=impuls(i)/w
        sdimp(i) = sqrt((sdimp(i)/w) - impuls(i)**2)
 2070 continue
c 
c-----------------------------------------------------------------------
c
      call barra('-',70,fout) 
      write(fout,2075)  w
 2075 format('>> sum of the weights:', f15.4)     
      rnum=2.0d0*numrep
      wm=w/rnum
      sdw=sqrt((sdw/rnum) - wm**2)
c
c-----------------------------------------------------------------------
c
      call barra('.',70,fout) 
      write(fout,*) '>> moments of the weights:'
      write(fout,*) '  -----------------------------------------' 
      write(fout,*) '           MEAN             STD'
      write(fout,*) '  -----------------------------------------' 
      write(fout,2080) wm,sdw
 2080 format(t8,2g15.4)
c
c-----------------------------------------------------------------------
c
      call barra('-',70,fout) 
      write(fout,*) '>> posterior moments of beta:'
      write(fout,*) '  -----------------------------------------' 
      write(fout,*) '        vname       MEAN      STD' 
      write(fout,*) '  -----------------------------------------' 
      do 2090 i=1+kdfix,np
 2090 write(fout,2091) name(i), betam(i), betasd(i)
 2091 format(6x,a8,2x,f8.4,2x,f8.4)
c
c-----------------------------------------------------------------------
c
      call barra('-',70,fout) 
      write(fout,*) '>> posterior moments of impulse responses:'
      call barra('-',65,fout) 
      write(fout,2099)
 2099 format(t2, 'I(i)',t15,'MEAN',t27,'STD',
     &       t36,'I(i)',t46,'MEAN',t60,'STD')
      call barra('-',65,fout) 
      do 2100 i=1,nimp/2
      write(fout,2110) i, impuls(i), sdimp(i),
     & nimp/2 + i, impuls(nimp/2 + i), sdimp(nimp/2 + i)
 2100 continue
 2110 format(t2,       'I(',i2,')',2x,f10.4,2x,f10.4,
     &       2x,'|',2x,'I(',i2,')',2x,f10.4,2x,f10.4)
c
c-----------------------------------------------------------------------
c
c compute means and sd's of 2 first predictive moments
c
      do 3000 i=1,12
      ypm(i)  = log(yT) + ypm(i)/w
 3000 ypsd(i) = ypsd(i)/w
c
      call barra('-',70,fout) 
      write(fout,*) '>> First 2 predictive moments of log(GNP_{T+i}):'
      write(fout,*)  '  ------------------------------------------'
      write(fout,*)  '  time               MEAN            STD'
      write(fout,*)  '  ------------------------------------------'
      do 3200 i=1,12
 3200 write(fout,3210) i, ypm(i), ypsd(i)
 3210 format(t4,'T+',i2,t18,f10.4,t33,f10.4)
      write(fout,*)  '  ------------------------------------------'
      write(fout,3215) log(yT)      
 3215 format('    (log(GNP_T) is: ',f10.4,' )')
      call barra('-',70,fout)
c-----------------------------------------------------------------------
      stop
      end
c
c***********************************************************************
c***********************************************************************
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c          end of main program --- subroutines and fcns follow
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c***********************************************************************
c***********************************************************************
c
c-----------------------------------------------------------------------
c  file handling
c-----------------------------------------------------------------------
c
       block data
        common/files/fdata,fpara,fout,ferr,fmma
        integer fdata,fpara,fout,ferr,fmma
        data fdata,fpara,fout,ferr,fmma/10,11,12,13,14/
       end
c
c-----------------------------------------------------------------------
c          subroutine: frac  (from sowell's code): changed
c-----------------------------------------------------------------------
c
c  subroutine to evaluate the log posterior for omega at a draw
c  for an arfima model. arfima(p,d,q)  d<.499
c
c-----------------------------------------------------------------------
      subroutine frac(theta,f,maxnpred,ierr)
c-----------------------------------------------------------------------
c
c   inputs to the routine are:
c
c    d = the differencing parameter   ( d < .5, but not 0)
c    p = the order of the autoregressive (ar)  polynomial (p<10)
c    q = the order of the moving average (ma) polynomial (p<10)
c    a = array of dimension 10 where the first p+1 elements are
c        the coefficients of the ar polynomial, (a(1) = 1)
c        a(1) + a(2)*z + a(3)*z**2 + a(4)*z**3 + ... + a(p+1)*z**p
c        the evaluation assumes the process is stationary hence
c        the roots of the ar polynomial must be outside the unit
c        circle.
c    b = array of dimension 10 where the first q+1 elements are
c        the coefficients of the ma polynomial, (b(1) = 1)
c        b(1) + b(2)*z + b(3)*z**2 + b(4)*z**3 + ... + b(q+1)*z**q
c
c   outputs from the routine are:
c
c      ierr = error flag if ierr = 0 after the routine then the
c           evaluation was successful otherwise a
c           problem was encountered.  the value of
c           ier will identify where the problem was.
c          1 = d > .499
c          3 = either p>9 or q>9
c          4 = ar polynomial has a root inside the unit circle
c             ( abs(rho(i))>.99 )
c          5 = inversion of the covariance matrix failed.
c    f = value of the log posterior for omega
c
c     the data is passed to the subroutine through a common statement
c     the number of observations is also passed through common statement
c    this program is written to not allow any more than 500 observations
c
      implicit double precision (a-h,o-z)
      complex zed(9),rho(9),hyp(1060,9),rho2p(9),sig1,ga(10)
      logical polish
      integer s,p,q
      integer fdata,fpara,fout,ferr,fmma
      dimension a(4),b(4),theta(7),c(530,9), rat(1040)
      common/data/x(500)
      common/parms/mynobs,p,q
      common/covar/sig(500),v2(500),chl(124750)
      common/files/fdata,fpara,fout,ferr,fmma
      common/hat/uhat,sigtmp,sse,dm1
c
c
c  intialize polynomials since autocov formula for d=0 uses
c  up to a(4) and b(4) regardless of p and q
c
      do 100 i=2,4
      a(i) = 0.0d0
      b(i) = 0.0d0
 100  continue
c
c  match the theta parameters to the parameters in the paper
c
      d = theta(1)
      a(1) = 1.0d0
      if (p.eq.0) goto 1010
      do 1000 i = 2,p+1
 1000 a(i) = theta(i)
 1010 b(1) = 1.0d0
      if (q.eq.0) goto 1030
      do 1020 i = 2,q+1
 1020 b(i) = theta(p+i)
 1030 continue
c
      ierr = 0
c
c   check to see if -5.0 < d < .499, var > .0001 and p,q<10.
c
      if ( (d.ge.(.499)) .or. (d.le.(-5.0)) ) then
       ierr = 1
       goto 9998
      endif
      if ( (p.ge.10) .or. (q.ge.10)) then
       ierr = 3
       goto 9998
      endif
c
c  jump to 1201 if d is 0 and use B-J formulae to get autocov's
c
      if (d.eq.0) goto 1201
c
c  if there are AR terms goto 1039
c
      if (p.ne.0) goto 1039
c      
c  code to compute cov when there are no AR terms
c
      cons = zgamma(1.0d0-2.0d0*d)/((zgamma(1.0d0-d))*(zgamma(d)))
      ichgvr = q + (mynobs+maxnpred) +1
      rat(q+ichgvr) = zgamma(d+q)/zgamma(1.0d0-d+q)
      do 1035 i=q+ichgvr-1,1,-1
 1035 rat(i) = rat(i+1)*(1.0d0-d+i-ichgvr)/(d+i-ichgvr)
      do 1038 s=1,(mynobs+maxnpred)
      sig(s) = 0.d0
      do 1037 n=1,q+1
      do 1036 m=1,q+1
 1036 sig(s) = sig(s)+b(n)*b(m)*rat(n-m-(s-1)+ichgvr)
 1037 continue
 1038 sig(s) = cons*sig(s)
      goto 1305
c
c  start here if there are AR terms
c
 1039  continue
c
c  code to calculate the covariance when there are AR terms
c
c    calculate the roots of the ar polynomial use numerical recipes routine.
c    note that the subroutine calculates the roots of the polynomial,
c    but the program uses the inverse of the roots.
c
c    routine requires input of complex coefficients, so create a new
c    variable, ga, which is a complex version of a
c
       do 1040 i=1,p+1
        ga(i)=cmplx(real(a(i)),0.)
 1040  continue
c
       polish=.true.
       call zroots(ga,p,rho,polish)
c
c  check to see that none of the roots of the ar polynomial are too 
c  close to the unit circle.	
c
      do 1050 i = 1,p
       if (rho(i).ne.0) then
       rho(i) = 1.0/rho(i)
       else
       ierr = 8
       goto 9998
      endif
       w = abs(rho(i))
       if (w.gt.(.995d0)) then
        ierr = 4
        goto 9998
       endif
 1050  continue
c
c    calculate the zeds
c     (calculate 1/zed then invert to get zed)
c
      do 1110 i = 1,p
      	zed(i) = rho(i)
      	do 1060 j = 1,p
 1060 	zed(i) = zed(i)*( 1.0d0 - rho(j)*rho(i) )
      	if (i.eq.1) goto 1080
      	do 1070 k = 1,i-1
 1070 	zed(i) = zed(i)*(rho(i) - rho(k))
 1080 	if (i.eq.p) goto 1100
      	do 1090 k = i+1,p
 1090 	zed(i) = zed(i)*(rho(i) - rho(k))
 1100 	zed(i) = 1.0d0/zed(i)
 1110 continue
c
c  to calculate the covariance terms we need to evaluate the
c  functions zed(i)*c(d,d,p+n-m-s,rho(i)).  in this program these values
c  will be denoted by the double array c(*,*).  to evaluate c(*,*) we
c  need to calculate hypergeometric functions f(d+h,1,1-d+h,rho(i)).
c  these values will be denoted by hyp(h-ichgvr,i). ichgvr is a
c  change of variables constant so that the indices on the arrays are
c  always positive.  ichgvr = p-q-mynobs.
c
c        first we need to evaluate all the hypergeometric
c        functions needed.  this will be done recursively
c        hence only one hypergeometric function needs to
c        be evaluated for each root.
c
      ichgvr = p-q-(mynobs+maxnpred)
      do 1130 i = 1,p
      call dhyper(d-ichgvr-1,1-d-ichgvr-1,rho(i),hyp((-1-2*ichgvr),i))
      j = (-1-2*ichgvr)
      do 1120 j = -2-2*ichgvr,1,-1
 1120 hyp(j,i) = dcmplx(1.0d0) + rho(i)*hyp(j+1,i)
     &      *dcmplx((d+j+ichgvr)/(1.0d0-d+j+ichgvr))
 1130 continue
c
c   now calculate the c(*,*) functions
c
      do 1140 i = 1,p
 1140 rho2p(i) = rho(i)**(2*p)
c
      coef=zgamma(1.0d0-2.0d0*d)/(zgamma(1.0d0-d)*zgamma(d))
c
      call ratio(d,ichgvr,r)
      r = coef*r
      do 1160 n = 1,2*q+(mynobs+maxnpred)
      	r = r*(d+ichgvr+n-1.0d0)/(1-d+ichgvr+n-1.0d0)
      	do 1150 i = 1,p
       	c(n,i) = zed(i)*r*( rho2p(i)*hyp(n,i) +
     &               hyp(-n-2*ichgvr,i) - 1.0d0 )
 1150 continue
 1160 continue
c
c    now calculate the autocovariances
c
      do 1200 s = 1,(mynobs+maxnpred)
        sig(s) = 0.0d0
        do 1190 i = 1,p
          sig1 = 0.0d0
          do 1180 n = 1,q+1
            do 1170 m = 1,q+1
 1170       sig1 = sig1 + b(n)*b(m)*c(p+n-m-(s-1)-ichgvr,i)
 1180     continue
          sig(s) = sig(s) + sig1
 1190   continue
 1200 continue
c 
 1201 continue
c
c if d = 0 
c get autocovariances using ARMA(3,3) formula
c
      if (d.eq.0) then
c 
       do 1210 i=1,p+1
        ga(i)=cmplx(real(a(i)),0.)
 1210  continue
c
       polish=.true.
       call zroots(ga,p,rho,polish)
c
c  check to see that none of the roots of the ar polynomial are too 
c  close to the unit circle.	
c
      do 1220 i = 1,p
       if (rho(i).ne.0) then
       rho(i) = 1.0/rho(i)
       else
       ierr = 8
       goto 9998
       endif
       w = abs(rho(i))
       if (w.gt.(.995d0)) then
        ierr = 4
        goto 9998
       endif
 1220  continue
c
c
      xx =  -1.0d0 - b(2)**2 - b(3)**2 + 2.0d0*a(2)**3*b(4) +    
     &   2.0d0*a(4)*b(4) - 2.0d0*a(4)**3*b(4) - b(4)**2 + 
     &   2.0d0*a(2)**2*(-b(3) + a(4)*b(4) - b(2)*b(4)) + 
     &   2.0d0*a(3)**2*(b(3) + a(4)*b(4) + b(2)*b(4)) + 
     &   a(4)**2*(1.0d0 + b(2)**2 + b(3)**2 + b(4)**2) + 
     &   a(3)*(-1.0d0 - b(2)**2 + 2*b(3) - b(3)**2 + 2.0d0*b(2)*b(4) - 
     &   b(4)**2 + 2.0d0*a(4)*(-b(2) - b(2)*b(3) + b(4) - b(3)*b(4))) + 
     &   a(2)*(b(2)*(2 + 2.0d0*b(3)) - 4.0d0*a(3)*b(4) - 
     &   2.0d0*a(3)**2*b(4) - 2.0d0*a(4)**2*b(4) + 2.0d0*b(3)*b(4) + 
     &   a(4)*(1.0d0 + b(2)**2 - 2.0d0*b(3) + b(3)**2 - 
     &   2.0d0*b(2)*b(4) + b(4)**2))
c            
      zz =  -1.0d0 + a(3)**3 + a(2)**3*a(4) + 2*a(4)**2 - a(4)**4 + 
     &   a(2)*a(4)*(1.0d0 - 4.0d0*a(3) - a(3)**2 - a(4)**2) + 
     &   a(3)*(-1.0d0 + a(4)**2) + a(3)**2*(1.0d0 + a(4)**2) + 
     &   a(2)**2*(1.0d0 - a(3) + a(4)**2)
c
      if (zz.eq.0) then
        write(ferr,*) 'denominator in sig(1) is 0'
        ierr = 8
        goto 9998
      endif       
      sig(1) = xx/zz
c
      xx = a(2) -a(3)*a(4) - b(2)- a(2)**2*b(2) +a(3)**2*b(2) +
     &   a(4)**2*b(2) +a(2)*b(2)**2- a(3)*a(4)*b(2)**2 + a(2)*b(3) +
     &   a(2)**3*b(3)-2.0d0*a(2)*a(3)*b(3)-a(2)*a(3)**2*b(3)+a(4)*b(3)+
     &   a(2)**2*a(4)*b(3)+ a(3)**2*a(4)*b(3) -  a(2)*a(4)**2*b(3) -
     &   a(4)**3*b(3)-b(2)*b(3)-a(2)**2*b(2)*b(3) +a(3)**2*b(2)*b(3) +
     &   a(4)**2*b(2)*b(3)+a(2)*b(3)**2-a(3)*a(4)*b(3)**2-a(2)**2*b(4)-
     &   a(2)**4*b(4) + a(3)*b(4) + 3.0d0*a(2)**2*a(3)*b(4) + 
     &   a(2)**2*a(3)**2*b(4) - a(3)**3*b(4) - 3.0d0*a(2)*a(4)*b(4) -
     &   a(2)**3*a(4)*b(4)-a(2)*a(3)**2*a(4)*b(4)+a(2)**2*a(4)**2*b(4)+
     &   a(3)*a(4)**2*b(4) +a(2)*a(4)**3*b(4) + a(2)*b(2)*b(4) +
     &   a(2)**3*b(2)*b(4)-2.0d0*a(2)*a(3)*b(2)*b(4)-
     &   a(2)*a(3)**2*b(2)*b(4) +
     &   a(4)*b(2)*b(4)+a(2)**2*a(4)*b(2)*b(4)+a(3)**2*a(4)*b(2)*b(4) -
     &   a(2)*a(4)**2*b(2)*b(4) - a(4)**3*b(2)*b(4) -  b(3)*b(4) -
     &   a(2)**2*b(3)*b(4) + a(3)**2*b(3)*b(4) + a(4)**2*b(3)*b(4) + 
     &   a(2)*b(4)**2 - a(3)*a(4)*b(4)**2 
      zz = 1.0d0 - a(2)**2 + a(3) +
     &   a(2)**2*a(3) - a(3)**2 - a(3)**3 - a(2)*a(4) -  a(2)**3*a(4) +
     &   4.0d0*a(2)*a(3)*a(4)+a(2)*a(3)**2*a(4)-2.0d0*a(4)**2-
     &   a(2)**2*a(4)**2 - 
     &   a(3)*a(4)**2 - a(3)**2*a(4)**2 + a(2)*a(4)**3 + a(4)**4
c
      if (zz.eq.0) then
        write(ferr,*) 'denominator in sig(2) is 0'
        ierr = 8
        goto 9998
      endif       
      sig(2) = -xx/zz
c
      sig(3) = b(3) - sig(2)*(a(2)+a(4)) - sig(1)*a(3) - a(2)*b(4) + 
     &         b(2)*b(4)
      sig(4) = b(4) - sig(3)*a(2) - sig(2)*a(3) - sig(1)*a(4) 
c
      do 1230 j=5,(mynobs+maxnpred)
 1230 sig(j) = - a(2)*sig(j-1) - a(3)*sig(j-2) - a(4)*sig(j-3)
      endif
c
c  the inverse and determinant of the covariance matrix is calculated
c  by first calculating the inverse and determinant of the correlation
c  matrix.
c  sig(s) contains sowell's gamma(s) ignoring the sigma-squared term
c  it corresponds to our v(i,j) elements
c
 1305  sigtmp = sig(1)
       if (sigtmp.ne.0) then
       do 1310 s = 1,(mynobs+maxnpred)
 1310  sig(s) = sig(s)/sigtmp
       else
       ierr = 7
       goto 9998
       endif
c
c  levin is a subroutine that calculates the log of the determinant and
c  the cholesky decomposition of the inverse of a toeplitz matrix.
c
       call levin(dlgdet,ier)
c
c  check to see if the levin subroutine was successful
c
      if (ier.ne.0) then
       ierr = 5
       goto 9998
      endif
c
c stuff before was for correlation matrix (call it rho)
c now convert back for determinant, for the rest be careful
c the fact that we are working with the correlation matrix instead of v11
c directly does not affect mu-hat, but we will have to rescale other 
c terms
c
      dlgdet=dlgdet+mynobs*dlog(sigtmp)
c
c calculate e_T' V^{-1}_{11} e_T
c
      dm1 = 1.0d0/v2(1) + ( (1.0d0 - chl(1))**2)/v2(2)
c
      do 1330 i = 3,mynobs
       z = 1.0d0
       do 1320 l = 1,i-1
 1320   z = z - chl( (i*(i-1)/2)+1-l)*1.0d0
 1330    dm1 = dm1 + z*z/v2(i)
c
c calculate e_T' V^{-1}_{11} w
c
        dm2= x(1)*1.0d0/v2(1) +
     &        (1.0d0 - chl(1))* (x(2) -chl(1)*x(1))/v2(2)
      do 1350 i = 3,mynobs
       z = 1.0d0
       z2=x(i)
       do 1340 l = 1,i-1
        z2=z2 - chl((i*(i-1)/2)+1-l)*x(l)
 1340   z = z - chl( (i*(i-1)/2)+1-l)*1.0d0
 1350   dm2 = dm2 + z*z2/v2(i)
c
c calculate mu-hat
c
       uhat=dm2/dm1
c
c sse
c
      sse=(x(1)-uhat)*(x(1)-uhat)/v2(1) +
     &    ((x(2) - uhat -chl(1)*(x(1)-uhat))**2)/v2(2)
      do 1370 i=3,mynobs
      z=x(i) - uhat
      do 1360 l=1,i-1
 1360 z=z - chl((i*(i-1)/2)+1-l)*(x(l)-uhat)
 1370 sse = sse + z*z/v2(i)
c
c now reconvert sse and dm1 to compensate for the fact that they were
c originally done in terms of the correlation matrix
c
      sse=sse/sigtmp
c
      dm1=dm1/sigtmp
c
c now add up all the individual terms of the log posterior
c
      f= -0.5d0*(dlgdet + log(dm1) + (mynobs-1)*log(sse))
c
c so far we have no prior, add log of the prior on here
c
      dlgpri=0.0d0
      f = f + dlgpri      
c
      goto 9999
 9998 continue
       if ((ierr.eq.4).or.(ierr.eq.8)) then
          write(ferr,*) 'ierr...',ierr, abs(rho(i))
         else  if (ierr.eq.7) then 
          write(ferr,*) 'ierr...',ierr, sigtmp
         else
          write(ferr,*) 'ierr...',ierr
       endif
      return
 9999 continue
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: frac
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------      
c  predictives' stuff
c-----------------------------------------------------------------------      
c      
c
      subroutine predict(npred,ypm,ypp)
c
      implicit double precision (a-h,o-z)
c      
      integer npred, mynobs,p,q
      double precision uhat, sigtmp, sse, dm1
      double precision z1, z2, zz1(500)  
      double precision ypm, ypp, prec0, prec1, prec2
c
      common/data/x(500)
      common/covar/sig(500),v2(500),chl(124750)
      common/parms/mynobs,p,q
      common/hat/uhat,sigtmp,sse,dm1
c
c
c  1. compute e'_n V_{21} = zz1
c
      do 2010 j=1,mynobs
      sum = 0.0d0
      do 2000 i=mynobs+1,mynobs+npred
 2000 sum = sum + sig(abs(i-j)+1)
 2010 zz1(j) = sum*sigtmp
c
c  2. compute [e'_n V_{21}] V^{-1}_{11} [w - muhat e_T]
c
      z = (zz1(1)*(x(1)-uhat))/v2(1)
      z = z+
     &   (zz1(2)-chl(1)*zz1(1))*(x(2)-uhat-chl(1)*(x(1)-uhat))/v2(2)
      do 2030 i=3,mynobs
      z1 = zz1(i)
      z2 = x(i)-uhat
      do 2020 j=1,i-1
      z1 = z1-chl((i*(i-1)/2)+1-j)*zz1(j)
 2020 z2 = z2-chl((i*(i-1)/2)+1-j)*(x(j)-uhat)
 2030 z = z + z1*z2/v2(i)
      z = z/sigtmp
c
c 3. compute location vector from last obs --- must add y(mynobs) later
c
      ypm = npred*uhat + z
c
c
c 4. compute e_n' V_{22} e_n = prec0
c    (It is just the sum of all elements in V_{22})
c
      sum = 0.0d0
      do 2040 j=1,npred
      do 2040 i=1,npred
 2040 sum = sum + sig(abs(i-j)+1)
      prec0 = sum*sigtmp
c
c 5. compute [e_n' V_{21}] V^{-1}_{11} [V_{12} e_n] = prec1
c
      z = (zz1(1)**2)/v2(1)
      z = z + ((zz1(2)-chl(1)*zz1(1))**2)/v2(2)
      do 2060 i=3,mynobs
      z1 = zz1(i)
      do 2050 j=1,i-1
 2050 z1 = z1-chl((i*(i-1)/2)+1-j)*zz1(j)
 2060 z = z + (z1**2)/v2(i)
      prec1 = z/sigtmp
c
c  make sure V_{22.1} is a s.p.d. matrix
c
c      if ((prec0-prec1).lt.0) write(*,*) 'V22.1 trouble!'
c         ierr = 9
c	 goto 3000
c      endif 
c
c 6. compute e_n' V_{21} V^{-1}_{11} e_T = prec2
c
      z = zz1(1)/v2(1)
      z = z + (zz1(2)-chl(1)*zz1(1))*(1.0d0-chl(1))/v2(2)
      do 2080 i=3,mynobs
      z1 = zz1(i)
      z2 = 1.0d0
      do 2070 j=1,i-1
      z1 = z1-chl((i*(i-1)/2)+1-j)*zz1(j)
 2070 z2 = z2-chl((i*(i-1)/2)+1-j)
 2080 z = z + z1*z2/v2(i)
      prec2 = z/sigtmp
c
c 7. put together pieces and get precision
c
      ypp=((mynobs-1)/sse)/(prec0-prec1+((npred-prec2)**2)/dm1)      
c
 3000 continue
      return
      end
c 
c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c          subroutine: levin (from sowell's code): unchanged
c-----------------------------------------------------------------------
c
c  this subroutine calculates "the log of the determinant"  and the 
c  cholesky decomposition of "the inverse" of a toeplitz matrix.  
c  
c  the log determinant is of the toeplitz matrix not the log 
c  determinant of the inverse of the toeplitz matrix
c
c-----------------------------------------------------------------------
c
      subroutine levin(dlgdet,err)
c
      implicit double precision (a-h,o-z)
      integer fdata,fpara,fout,ferr,fmma
      integer p,q,err
c
      common/covar/sig(500),v2(500),a(124750)
      common/parms/mynobs,p,q
      common/files/fdata,fpara,fout,ferr,fmma
c
      err = 0                  
      v2(1) = sig(1)
      if (v2(1).le.0.0) then
         write(ferr,*) 'error in levin, sig(1) is not positive'
         err = 1
         goto 50
         endif
      a(1) = sig(2)/sig(1)
      v2(2) = v2(1)*(1.0d0 - a(1)**2)
      if (v2(2).le.0.0) then
         write(ferr,*) 'error in levin, v2(2) is not positive'
c	 write(ferr,*) sig(1), sig(2)
         err = 2
         goto 50
         endif
      do 30 i = 2,mynobs-1
         temp = 0.0d0
         do 10 k = 1,i-1
10          temp = temp + a(((i-1)*(i-2)/2)+k)*sig(i-k+1)
         a((i*(i+1)/2)) = (sig(i+1) - temp)/v2(i)
         do 20 j = 1,i-1
20          a((i*(i-1)/2)+j) = a(((i-1)*(i-2)/2)+j) -
     &            a((i*(i+1)/2))*a((i*(i-1)/2)+1-j)
         v2(i+1) = v2(i)*(1.0d0 - a((i*(i+1)/2))**2)
      if (v2(i+1).le.0.0) then
         write(ferr,*) 'error in levin, cov not p.d.', v2(i+1)
         err = 3
         goto 50
         endif
30    continue
      dlgdet = 0.0d0
      do 40 i = 1,mynobs
40        dlgdet = dlgdet + log(v2(i))
50    return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: levin
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: dhyper (from sowell's code): unchanged
c-----------------------------------------------------------------------
c
c  this is a subroutine to evaluate a hypergeometric function of
c  the form  f(a,1;c,z).  the approach is to simply sum the series
c  until the terms are below 1.0d-14 or until the first 10,000 terms
c  have been summed
c
c  note: the hypergeometric function converges very slowly when z is
c  near the unit circle.
c
      subroutine dhyper(a,c,z,hyp)
      implicit double precision (a-h,o-z)
      integer fdata,fpara,fout,ferr,fmma
      complex z,hyp,inc
      common/files/fdata,fpara,fout,ferr,fmma

c
      acc = 1.0d-14
      hyp = (1.0d0,0.0d0)
      inc = (1.0d0,0.0d0)
      do 10 m = 1,10000
       inc = inc*z*(a+m-1.0d0)/(c+m-1.0d0)
       hyp = hyp + inc
       w = abs(inc)
 10     if (w.lt.acc) goto 20
       write(ferr,*) 'all through the hyper loop !!!  10,000 iterations'
       write(ferr,*) 'a,c,z,hyp are :',a,c,z,hyp
 20    return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: dhyper
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: zgamma (from sowell's code): changed 
c-----------------------------------------------------------------------
c
c this is my gamma function which takes the basic function (dgamma)
c and extends it to allow for negative arguments
c note that it does not handle properly arguments of zero or negative
c integers (these imply infinite values for the gamma function)
c but this is no problem unless d=0 (and we rule out this case)
c
      double precision function zgamma(xx)
c
      implicit double precision (a-h,o-z)
      integer lint,j
c
      if (xx.gt.0) then
       zgamma=gamma(xx)
       goto 5
      endif
      lint=0
      lint=int(xx)
      lint=-lint
      dec=-(xx+lint)
      xx1=1-dec
      zgamma=gamma(xx1)
      tmp1=xx
      do 12 j=1,lint
       tmp1=tmp1*(xx+j)
 12    continue
      zgamma=zgamma/tmp1
 5    continue
      return
      end
c
c-----------------------------------------------------------------------
c          end of function: zgamma
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: ratio (from sowell's code): unchanged YET
c-----------------------------------------------------------------------
c
c   this subroutine calculates the ratio of two gamma functions of the 
c   form gamma(d+k)/gamma(1-d+k).  the problem is that gamma(x) is 
c   machine infinite on systems typically at x=35 or x=56 or x=157.
c   note that
c   
c     gamma(d+k)       gamma(d-k)
c   -------------- = --------------
c    gamma(1-d+k)     gamma(1-d-k)
c
c-----------------------------------------------------------------------
c
      subroutine ratio(d,k,r)
c
      implicit double precision (a-h,o-z)
      r = zgamma(d)/zgamma(1-d)
      j = abs(k)
      do 10 i = 1,j
 10     r = r*(d+i-1.0d0)/(1.0d0-d+i-1.0d0)
      return
      end
c-----------------------------------------------------------------------
c          end of subroutine: ratio
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: zimp
c-----------------------------------------------------------------------
c
c calculate the impulse responses up to order nimp
c since we are basing things on sowell's code, 
c note that signs are reversed
c of the beta's and phi's so we must be careful
c
      subroutine zimp(beta,nimp,din)
c
      implicit double precision (a-h,o-z)
      integer p,q
      dimension beta(7),din(41),f(101),thma(11),djj(101)
      common/parms/mynobs,p,q
c
c first get the f's (we need f(h) for h=1,..,nimp+1)
c
      f(1)=1.0d0
      do 217 i=2,nimp+1
        f(i)=0.0d0
        do 216 j=1,p
          i1=i-j
          if (i1.ge.1) then
            f(i)= f(i) - beta(j+1)*f(i1)
          endif
 216    continue
 217  continue
c 
c now calculate the js
c
      djj(1)=1.0d0
      thma(1)=1.0d0
      do 218 i=1,q
 218    thma(i+1)=beta(p+1+i)
      do 220 i=2,nimp+1
        djj(i)=0.0d0
        do 219 j=1,q+1
          i1=i-j+1
          if (i1.gt.0) then
             djj(i)=djj(i) + thma(j)*f(i1)
          endif
 219    continue
 220  continue
c
c calculate the impulse responses
c
      do 221 i=1,nimp
        din(i) = djj(i+1)
        ci=1.0d0
        do 222 j=1,i
          realj=j
          ci=ci*(realj+beta(1))/realj
          din(i)=din(i) + ci*djj(i+1-j)
 222    continue
 221  continue
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: zimp
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: laguer
c-----------------------------------------------------------------------
c
c laguer is called by zroots to take the roots of a polynomial
c drawn from numerical recipes p. 264-265
c
      subroutine laguer(a,m,x,eps,polish)
c
c the root finding subroutine taken from numerical recipes
c
      complex a(m+1),x,dx,x1,b,d,f,g,h,sq,gp,gm,g2,zero
      integer iter,j
      logical polish
      parameter (zero=(0.,0.),epss=6.e-8,maxit=100)
      dxold=cabs(x)
      do 12 iter=1,maxit
	   b=a(m+1)
	   err=cabs(b)
	   d=zero
	   f=zero
	   abx=cabs(x)
	   do 11 j=m,1,-1
	    f=x*f+d
	    d=x*d+b
	    b=x*b+a(j)
	    err=cabs(b)+abx*err
 11	   continue
	   err=epss*err
	   if(cabs(b).le.err) then
	    dx=zero
	    return
	   else
	    g=d/b
	    g2=g*g
	    h=g2-2.*f/b
	    sq=csqrt((m-1)*(m*h-g2))
	    gp=g+sq
	    gm=g-sq
	    if(cabs(gp).lt.cabs(gm)) gp=gm
	    dx=m/gp
	   endif
	   x1=x-dx
	   if(x.eq.x1) return
	   x=x1
	   cdx=cabs(dx)
	   if(.not.polish)then
	    if(cdx.le.eps*cabs(x)) return
	   endif
 12   continue
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: laguer
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: zroots
c-----------------------------------------------------------------------
c
c finds roots of polynomial, drawn from numerical recipes p. 265-266
c
      subroutine zroots(a,m,roots,polish)
c
      parameter (eps=1.e-6,maxm=101)
      integer j,jj,i
      complex a(m+1),roots(m),ad(maxm),x,b,c
      logical polish
      do 11 j=1,m+1
	   ad(j)=a(j)
 11    continue
      do 13 j=m,1,-1
	   x=cmplx(0.,0.)
	   call laguer(ad,j,x,eps,.false.)
       if (abs(aimag(x)).le.2.*eps**2*abs(real(x))) x=cmplx(real(x),0.)
	   roots(j)=x
	   b=ad(j+1)
	   do 12 jj=j,1,-1
	    c=ad(jj)
	    ad(jj)=b
	    b=x*b+c
 12	   continue
 13   continue
      if (polish) then
	   do 14 j=1,m
	    call laguer(a,m,roots(j),eps,.true.)
 14	   continue
      endif
      do 16 j=2,m
	   x=roots(j)
	   do 15 i=j-1,1,-1
	    if(real(roots(i)).le.real(x)) go to 10
	    roots(i+1)=roots(i)
 15	   continue
       i=0
 10    roots(i+1)=x
 16   continue
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: zroots
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: lngamma
c-----------------------------------------------------------------------
c
c This function computes ln(gamma(x)) for x > 0.
c Taken from Pascal Numerical Recipes p 177
c
c Modified to improve accuracy when x<1 using:
c
c                     pi*e
c Gamma(1-e) = ---------------------
c               Gamma(1+e)*sin(pi*e)
c
      double precision function gammln(xx)
c
      double precision xx
      double precision cof(6),stp,half,one,fpf,x,y,tmp,ser,g,dpi
c
      data cof/76.18009173d0,-86.50532033d0,24.01409822d0, 
     &         -1.231739516d0,0.120858003d-2,-.536382d-5/
      data stp/2.50662827465d0/
      data dpi/3.14159265358979d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
c                 
      if (xx.lt.1) then
      x = one - xx
      tmp =  x + fpf
      tmp = (x + half)*log(tmp)-tmp
      ser = one
      do 1000 j=1,6
         x=x+one
         ser = ser+cof(j)/x
 1000 continue
      g = tmp + log(stp*ser)
      y = dpi*(one-xx)
      gammln = log(y/sin(y)) - g
c
      else
c
      x = xx - one
      tmp =  x + fpf
      tmp = (x + half)*log(tmp)-tmp
      ser = one
      do 1100 j=1,6
         x=x+one
         ser = ser+cof(j)/x
 1100 continue
      gammln = tmp + log(stp*ser)
      endif
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of function gammln
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: gamma
c-----------------------------------------------------------------------
c
c This function computes gamma(x) for x>0.
c
       double precision function gamma(x)
       double precision x,gammln
c
       gamma = exp(gammln(x))
c
       return
       end
c
c-----------------------------------------------------------------------
c          end of function gamma
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c     subroutine barra
c-----------------------------------------------------------------------
c
c
      subroutine barra(car,n,fout)
c
c-----------------------------------------------------------------------
c
      integer i, n, fout
      character*1 car
      character*1 bar(80)
      do 1000 i=1,n
 1000 bar(i) = car
      write(fout,*) (bar(i),i=1,n)
      return
      end
c
c-----------------------------------------------------------------------
c   end of subroutine barra
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          function: ran2
c-----------------------------------------------------------------------
c
      double precision function ran2(idum)
c
c   this function contains a portable random number generator 
c   uniform [0,1]
c   (numerical recipes p. 273)
c
      implicit double precision (a-h,o-z)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      double precision AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1.d0/IM1,IMM1=IM1-1,
     &IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     &NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2d-7,RNMX=1.0d0-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END      
c
c-----------------------------------------------------------------------
c          end of function: ran2
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c    fcn u draws uniformly in [lower,upper]
c-----------------------------------------------------------------------
c
      double precision function u(lower, upper)
c
      implicit double precision (a-h,o-z)
      double precision lower, upper
      common/seed/idum
c      
      u = (ran2(idum)*(upper-lower))+lower
c      
      return
      end
c
c-----------------------------------------------------------------------
c          subroutine: u1 
c-----------------------------------------------------------------------
c
      subroutine u1(draw1, draw2, ndraws)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision draw1(7), draw2(7)
      common/seed/idum
c
c draw beta1
c
      draw1(2) = u(-0.99d0,0.99d0)
      draw2(2) = -draw1(2)
      ndraws = ndraws+1
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u1 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u11 
c-----------------------------------------------------------------------
c
      subroutine u11(draw1, draw2, ndraws)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision draw1(7), draw2(7)
      common/seed/idum
c
c draw beta1
c
      draw1(2) = u(-0.99d0,0.99d0)
      draw2(2) = -draw1(2)
c
c draw beta2
c
      draw1(3) = u(-0.99d0,0.99d0)
      draw2(3) = -draw1(3)
c
      ndraws = ndraws+1
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u11 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c     subroutine drawtrian
c     draws uniformly in the stat/inv triangle given by
c         t1 < 1+t2 
c     and 
c         t1 > -1-t2
c-----------------------------------------------------------------------
c
      subroutine drawtrian(draw1, draw2, ndraws)
c
      implicit double precision (a-h, o-z)
      integer ndraws
      double precision draw1(2), draw2(2), t1, t2
      common/seed/idum
c
 1000 t1 = u(-2.0d0,2.0d0)
      t2 = u(-1.0d0,1.0d0)
      ndraws = ndraws + 1
      if ((t1.ge.(1+t2)).or.(t1.le.(-1-t2))) goto 1000
c
      draw1(1) = t1
      draw1(2) = t2
c
      draw2(1) = -t1
      draw2(2) =  t2
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: drawtrian
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c          subroutine: u2 
c-----------------------------------------------------------------------
c
      subroutine u2(draw1, draw2, ndraws)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(2), d2(2), draw1(7), draw2(7)
      common/seed/idum
c
c draw beta2 and beta3
c
      call drawtrian(d1,d2,ndraws)
c
      draw1(2) = d1(1)
      draw1(3) = d1(2)
c
      draw2(2) = d2(1)
      draw2(3) = d2(2)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u2 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u12 
c-----------------------------------------------------------------------
c
      subroutine u12(draw1, draw2, ndraws)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(2), d2(2), draw1(7), draw2(7)
      common/seed/idum
c
c draw beta2
c
      draw1(2) = u(-0.99d0,0.99d0)
      draw2(2) = -draw1(2)
c
c draw beta3 and beta4
c
      call drawtrian(d1,d2,ndraws)
c
      draw1(3) = d1(1)
      draw1(4) = d1(2)
c
      draw2(3) = d2(1)
      draw2(4) = d2(2)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u12 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u21 
c-----------------------------------------------------------------------
c
      subroutine u21(draw1, draw2, ndraws)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(2), d2(2), draw1(7), draw2(7)
      common/seed/idum
c
c draw beta2 and beta3
c
      call drawtrian(d1,d2,ndraws)
c
      draw1(2) = d1(1)
      draw1(3) = d1(2)
c
      draw2(2) = d2(1)
      draw2(3) = d2(2)
c
c draw beta4
c
      draw1(4) = u(-0.99d0,0.99d0)
      draw2(4) = -draw1(4)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u21 
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c          subroutine: u22 
c-----------------------------------------------------------------------
c
      subroutine u22(draw1, draw2, ndraws)
c
      implicit double precision (a-h, o-z)
      integer ndraws
      double precision draw1(7), draw2(7), d1(2), d2(2)
      common/seed/idum
c
c  AR---beta2 and beta3
c
      call drawtrian(d1, d2, ndraws)
      draw1(2) = d1(1)
      draw2(2) = d2(1)
      draw1(3) = d1(2)
      draw2(3) = d2(2)
c
c  MA---beta4 and beta5
c
      call drawtrian(d1, d2, ndraws)
c
c  since ndraws was incremented above already
c
      ndraws = ndraws - 1
c
      draw1(4) = d1(1)
      draw2(4) = d2(1)
      draw1(5) = d1(2)
      draw2(5) = d2(2)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: u22 
c-----------------------------------------------------------------------
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c          subroutine: drawtres
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c See: J. Monaham (1984): "A Note on Enforcing Stationarity in 
c                          ARMA Models", Biometrika, 71, 403-404
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine drawtres(draw1, draw2, zlit1, zlit2)
c
      implicit double precision (a-h, o-z)
      dimension r1(3), r2(3), draw1(3), draw2(3)
      common/seed/idum
c
c draw r's in [-1,1]
c
  10  continue
      r1(1) = u( -1.0d0, 1.0d0)
      r2(1) = -r1(1)
      r1(2) = u( -1.0d0, 1.0d0)
      r2(2) = -r1(2)
      r1(3) = u( -1.0d0, 1.0d0)
      r2(3) = -r1(3)
c
c transform r's into theta's and phi's
c
      draw1(1) = r1(1)+r1(2)*r1(3)+r1(1)*r1(2)
      draw1(2) = r1(2)+r1(3)*(r1(1)+r1(1)*r1(2))
      draw1(3) = r1(3)
c
c importance fcn stuff
c
      a = 1.0d0-r1(3)**2
      b = 1.0d0+r1(2)
      if ((a.gt.0.0d0).and.(b.gt.0.0d0)) then
       zlit1 = log(a) + log(b)
      else
       goto 10
      endif
c
c Same for the antithetic draw
c
c transform r's into theta's and phi's
c
      draw2(1) = r2(1)+r2(2)*r2(3)+r2(1)*r2(2)
      draw2(2) = r2(2)+r2(3)*(r2(1)+r2(1)*r2(2))
      draw2(3) = r2(3)
c
c importance fcn stuff
c
      a = 1.0d0-r2(3)**2
      b = 1.0d0+r2(2)
c
c  should never happen, but just in case:
c
      if ((a.gt.0.0d0).and.(b.gt.0.0d0)) then
       zlit2 = log(a) + log(b)
      else
       goto 10
      endif
c
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c          end of  subroutine: drawtres
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c-----------------------------------------------------------------------
c          subroutine: u3 
c-----------------------------------------------------------------------
c
      subroutine u3(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(3), d2(3), draw1(7), draw2(7)
      common/seed/idum
c
c draw beta2, beta3 and beta4
c
      call drawtres(d1,d2, zlit1, zlit2)
      ndraws = ndraws + 1
c
      draw1(2) = d1(1)
      draw1(3) = d1(2)
      draw1(4) = d1(3)
c
      draw2(2) = d2(1)
      draw2(3) = d2(2)
      draw2(4) = d2(3)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u3 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u13 
c-----------------------------------------------------------------------
c
      subroutine u13(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(3), d2(3), draw1(7), draw2(7)
      common/seed/idum
c
      ndraws = ndraws+1
c
c draw beta2
c
      draw1(2) = u(-0.99d0,0.99d0)
      draw2(2) = -draw1(2)
c
c draw beta3, beta4, beta5
c
      call drawtres(d1,d2, zlit1, zlit2)
      ndraws = ndraws + 1
c
      draw1(3) = d1(1)
      draw1(4) = d1(2)
      draw1(5) = d1(3)
c
      draw2(3) = d2(1)
      draw2(4) = d2(2)
      draw2(5) = d2(3)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u13 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u32 
c-----------------------------------------------------------------------
c
      subroutine u32(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(2), d2(2), draw1(7), draw2(7)
      double precision dd1(3), dd2(3)
      common/seed/idum
c
c draw beta2, beta3 and beta4
c
      call drawtres(dd1,dd2, zlit1, zlit2)
c
      draw1(2) = dd1(1)
      draw1(3) = dd1(2)
      draw1(4) = dd1(3)
      draw2(2) = dd2(1)
      draw2(3) = dd2(2)
      draw2(4) = dd2(3)
c
c draw beta5 and beta6
c
      call drawtrian(d1,d2,ndraws)
c
      draw1(5) = d1(1)
      draw1(6) = d1(2)
c
      draw2(5) = d2(1)
      draw2(6) = d2(2)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u32 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u31 
c-----------------------------------------------------------------------
c
      subroutine u31(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(3), d2(3), draw1(7), draw2(7)
      common/seed/idum
c
c draw beta2, beta3, beta4
c
      call drawtres(d1,d2, zlit1, zlit2)
      ndraws = ndraws + 1
c
      draw1(2) = d1(1)
      draw1(3) = d1(2)
      draw1(4) = d1(3)
c
      draw2(2) = d2(1)
      draw2(3) = d2(2)
      draw2(4) = d2(3)
c
c draw beta5
c
      draw1(5) = u(-0.99d0,0.99d0)
      draw2(5) = -draw1(4)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u31 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: u23 
c-----------------------------------------------------------------------
c
      subroutine u23(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer ndraws
      double precision d1(2), d2(2), draw1(7), draw2(7)
      double precision dd1(3), dd2(3)
      common/seed/idum
c
c draw beta2 and beta3
c
      call drawtrian(d1,d2,ndraws)
c
      draw1(2) = d1(1)
      draw1(3) = d1(2)
c
      draw2(2) = d2(1)
      draw2(3) = d2(2)
c
c draw beta4, beta5 and beta6
c
      call drawtres(dd1,dd2, zlit1, zlit2)
c
      draw1(4) = dd1(1)
      draw1(5) = dd1(2)
      draw1(6) = dd1(3)
c
      draw2(4) = dd2(1)
      draw2(5) = dd2(2)
      draw2(6) = dd2(3)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of  subroutine: u23 
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c          subroutine: u33 
c-----------------------------------------------------------------------
c
      subroutine u33(draw1, draw2, ndraws, zlit1, zlit2)
c
      implicit double precision (a-h, o-z)
      integer ndraws
      double precision draw1(7), draw2(7), dd1(3), dd2(3)
      common/seed/idum
c
      ndraws = ndraws + 1
c
c  AR---beta2 , beta3 and beta4
c
      call drawtres(dd1,dd2, zlit1, zlit2)
c
      draw1(2) = dd1(1)
      draw1(3) = dd1(2)
      draw1(4) = dd1(3)
c
      draw2(2) = dd2(1)
      draw2(3) = dd2(2)
      draw2(4) = dd2(3)
c
c  MA---beta5, beta6 and beta7
c
      call drawtres(dd1,dd2, zlit1, zlit2)
c
      draw1(5) = dd1(1)
      draw1(6) = dd1(2)
      draw1(7) = dd1(3)
c
      draw2(5) = dd2(1)
      draw2(6) = dd2(2)
      draw2(7) = dd2(3)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: u33 
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c          subroutine: draw
c-----------------------------------------------------------------------
c
      subroutine draw(draw1, draw2, dfixed, dvalue, p, q, ndraws,
     &                zlit1, zlit2)
c
      implicit double precision (a-h,o-z)
      integer p, q, ndraws
      double precision draw1(7), draw2(7), dvalue
      logical dfixed
      common/seed/idum
c
      zlit1 = 0.0d0
      zlit2 = 0.0d0
c
c
c beta1--- draw d in [-1.,.5]
c
      if (dfixed) then
        draw1(1) = dvalue
        draw2(1) = dvalue
      else
        draw1(1) = u(-1.0d0,0.5d0)
        draw2(1) = -0.25d0-(draw1(1)-(-0.25d0))
      endif
c
c beta2...7
c
      if (p.eq.0) then
        if (q.eq.0) then
            ndraws=ndraws+1
	    return
	else if (q.eq.1) then
	    call u1(draw1, draw2, ndraws)
	else if (q.eq.2) then
	    call u2(draw1, draw2, ndraws)
	else if (q.eq.3) then
	    call u3(draw1, draw2, ndraws, zlit1, zlit2)
	endif
      else if (p.eq.1) then
        if (q.eq.0) then
	    call u1(draw1, draw2, ndraws)
	else if (q.eq.1) then
	    call u11(draw1, draw2, ndraws)
	else if (q.eq.2) then
	    call u12(draw1, draw2, ndraws)
	else if (q.eq.3) then
	    call u13(draw1, draw2, ndraws, zlit1, zlit2)
	endif
      else if (p.eq.2) then	
        if (q.eq.0) then
	    call u2(draw1, draw2, ndraws)
	else if (q.eq.1) then
	    call u21(draw1, draw2, ndraws)
	else if (q.eq.2) then
	    call u22(draw1, draw2, ndraws)
	else if (q.eq.3) then
	    call u23(draw1, draw2, ndraws, zlit1, zlit2)
	endif
      else if (p.eq.3) then	
        if (q.eq.0) then
	    call u3(draw1, draw2, ndraws, zlit1, zlit2)
	else if (q.eq.1) then
	    call u31(draw1, draw2, ndraws, zlit1, zlit2)
	else if (q.eq.2) then
	    call u32(draw1, draw2, ndraws, zlit1, zlit2)
	else if (q.eq.3) then
	    call u33(draw1, draw2, ndraws, zlit1, zlit2)
	endif
      endif      
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of subroutine: draw
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c          function: ratpol (evaluates MA/AR ev at 1)
c-----------------------------------------------------------------------
c
      double precision function ratpol(beta, p, q)
c
      implicit double precision (a-h,o-z)
      integer p, q
      double precision sabove, sbelow, beta(7)
c      
      sabove = 1.0d0
      sbelow = 1.0d0
c  
      if (p.ne.0) then
      do 1000 i=1,p
 1000 sbelow = sbelow + beta(i+1)
      endif
c
      if (q.ne.0) then
      do 2000 j=1,q
 2000 sabove = sabove + beta(j+1+p)
      endif
c 
      ratpol = sabove/sbelow
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of function: ratpol (evaluates MA/AR ev at 1)
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c     subroutine acc
c-----------------------------------------------------------------------
c
      subroutine acc(datum,weight,datmin,datmax,nbins,cdf,b,a,su,sq)
c      
c-----------------------------------------------------------------------
c
      implicit double precision(a-h,o-z)
c      
      double precision datum, weight, datmin, datmax, b, a, su, sq
      integer i, nbins
      double precision cdf(nbins)
c
      double precision x, step
c
      su = su + weight*datum
      sq = sq + weight*datum**2
c
      if (datum.lt.datmin) b = b + weight
      if (datum.ge.datmax) a = a + weight
c
      x = datmin
      step = (datmax-datmin)/nbins
c
      do 1000 i=1,nbins
        x = x + step
        if (datum.le.x) cdf(i)=cdf(i)+weight
 1000 continue
      return
      end 
c
c-----------------------------------------------------------------------
c     end of subroutine acc
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c
      subroutine hist(name,letra,cdf,pdf,datmin,datmax,nbins,b,a,su,sq)
c
c-----------------------------------------------------------------------
c
      implicit double precision(a-h,o-z)
c 
      character*(*) name
      character*1 letra
      integer nbins, fdata,fpara,fout,ferr,fmma
      double precision datmin, datmax, b, a, su, sq, prev
      dimension cdf(nbins), pdf(nbins)
c
      common/files/fdata,fpara,fout,ferr,fmma
c      
      double precision max, sum, bound, step, midpoint
      integer i
c
      character*50 height
c 
      height='XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
c      
      step = (datmax-datmin)/nbins
      pdf(1) = cdf(1)
      sum = a + pdf(1)
      do 1000 i=2,nbins
      pdf(i) = cdf(i)-cdf(i-1)
      sum = sum + pdf(i)
 1000 continue
c
      su = su/sum
      sq = sqrt((sq/sum)-su**2)
c
c normalize
c
      max = 0.0d0
      do 2000 i=1,nbins
      pdf(i) = pdf(i)/sum
      if (pdf(i).gt.max) max=pdf(i)
 2000 continue 
      bound = datmin
      write(fout,*)
      write(fout,2060)
      write(fout,2050) name
 2050 format(t15,'%% ',t30, a10, t65,'%%')
       write(fout,2051)
 2051 format(t15, '%% ', t65, '%%')
      write(fout,2052) su,sq
 2052 format(t15,'%%   Mean...',g13.4,4x,'STD...',g10.4,t65, '%%')
      write(fout,2060)
 2060 format(t15,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
      write(fout,*)
c
      write(fout,2100) '   UBound    PDF   CDF(%)'
      write(fout,2100) '+---------------------------'
 2100 format(t51,a26)
      prev = 0.0d0
      do 3100 i=1,nbins
      bound = bound + step
      if (((cdf(i)/sum).ge.(0.001d0)).and.(prev.lt.(0.999d0))) then
      write(fout,3099) height(1:nint((pdf(i)/max)*50)), bound,
     & pdf(i), (cdf(i)/sum)*100.0d0
 3099 format(a50, '|',f9.3, 1x, f6.4, 1x, f8.4)
      endif
      prev = cdf(i)/sum
 3100 continue
c
c mass to the left and to the right of range bounds
c
      write(fout,*)
      write(fout,4000) datmin, (b/sum)*100.0d0
 4000 format(t15,'Mass to the left of  ',f16.4,' is.....',f8.4,' %')
      write(fout,4010) datmax, (a/sum)*100.0d0
 4010 format(t15,'Mass to the right of ',f16.4,' is.....',f8.4,' %')
c
c print out a table for Mathematica
c
      write(fmma,*) name//letra,' ={' 
      midpoint = datmin + step/2.0d0
      do 5100 i=1,nbins-1
      write(fmma,4099) midpoint,pdf(i)
 4099 format('{',f8.4,','f14.10,'},')
 5100 midpoint = midpoint + step
      write(fmma,4100) midpoint,pdf(nbins)
 4100 format('{',f8.4,','f14.10,'}};')
c      
      return
      end
c
c-----------------------------------------------------------------------
c     end of sub hist
c-----------------------------------------------------------------------
c                real FUNCTION secnds(t)
c                secnds = (LONG(362)/60.0)-t
c                END
c
c-----------------------------------------------------------------------
c          function: studt
c-----------------------------------------------------------------------
c
      double precision function studt(x,df,mean,prec)
c
      integer df
      double precision x,mean,prec
c
      studt=sqrt(abs(prec))*(1+(prec/df)*(x-mean)**2)**(-(df+1)*0.5d0)
c
      return
      end
c
c-----------------------------------------------------------------------
c          end of function: studt
c-----------------------------------------------------------------------
c

	
    
	

     
