c
c
c*************************************************************************
c
c
      subroutine STAT ( m, n, nkg2, kg, y, weight, c , alf, a,
c
     *                  iv, t, eta, r, 
c
     *                  tlabel, tlabld, ylabel, ylabld )
c
c------------------------------------------------------------------------
c
c   provide the indicated information in response to
c   the rs command: 
c  (1)  standard deviation of the parameters
c  (2)  t-ratio
c  (3)  residual sum of squares
c  (4)  residual standard error
c  (5)  coefficient of determination
c  (6)  standardized residuals
c  (7)  plot of standardized residuals vs indendent
c       variables
c  (8)  plot of standardized residuals vs predicted
c       values
c  (9)  normal probability plot.
c
c-------------------------------------------------------------------------------
c
      PARAMETER ( mmax = 3000, maxn = 100, maxkg = 50, miv = 6 )
c
      character*1 ians
      character*47  tlabel(miv), ylabel
      character*48  tlabld(miv), ylabld
c
      real  y(m), c(maxn), alf(maxkg), a(mmax,nkg2)
c
      real  t(mmax,iv), weight(m), eta(m), r(m), temp(mmax)
c
      real  bm(mmax,maxn+maxkg), xx1(mmax,maxn+maxkg), ph(maxn+1)
c
      real  syh(mmax), sr(mmax), srv(mmax) 
c
      logical  ip, ts, idp, idcp, itp
c
      dimension  ior(mmax)
c
c
      common /blokip/ ip
      common /blkph/ ph
      common /blkbm/ bm
      common /blkts/ ts
      common /blkitp/ itp
      common /blkidp/ idp
      common /blidcp/ idcp
c
c-------------------------------------------------------------------------------
c
c--  Initialize values
c
      npkg=n+kg
      rms=a(2,nkg2) 
      pi = acos( -1.0 )
c
c--  Print parameters, their std dev & t-ratios
c
      print 80
      if (ip) write (10,80)
80    format ('1'/' '/4x,
     1            'parameter',6x,'s.d. of parm.',5x,'t-ratio'/' ')
      if (n.eq.0) go to 110
      do 100 i=1,n
         ai=abs(a(i,i))
         bb=sqrt(ai)
         tr=c(i)/bb 
         print 90, c(i),bb,tr 
         if (ip) write (10,90) c(i),bb,tr
90       format (' ',1p3e16.6)
100   continue
110   if (kg.eq.0) go to 130
      do 120 i=1,kg 
         j=n+i
         aj=abs(a(j,j))
         bb=sqrt(aj)
         tr=alf(i)/bb
         print 90, alf(i),bb,tr
120      if (ip) write (10,90) alf(i),bb,tr
c
c--  Print residual sum of squares, res mean square, res std error, and
c          coefficient of determination.
c
 130  continue
c
      rss=(a(1,nkg2))**2
      rsd=sqrt(abs(rms))
      print 140, rss
      if (ip) write (10,140) rss
140   format (' '/10x, ' wtd residual sum of squares:  ',1pe16.6)
      print 150, rms
      if (ip) write (10,150) rms
150   format (' '/' ',10x,'wtd residual mean square:  ',3x,1pe16.6)
      print 160, rsd
      if (ip) write (10,160) rsd
160   format (' '/' ',10x,'wtd residual standard error:  ',1pe16.6/' ')
      yb=0.
      sws=0.
      do 170 i=1,m
         sws=sws+weight(i)
170      yb=yb+y(i)*weight(i) 
      yb=yb/sws
      ctss=0.
      do 180 i=1,m
180      ctss=ctss+weight(i)*(y(i)-yb)**2
      r2=1.-rss/ctss
      print 190, r2 
      if (ip) write (10,190) r2
190   format (' coefficient of determination (r-square): ',1pe16.6/' ')
c
c-------------------------------------------------------------------------------
c
c--  Calculate standardized residuals
c
200   if (n.eq.0) go to 220
      do 210 i=1,m
         ii=i
c
c       -----------------------------------------------------
         call CALC ( ii, kg, n, mmax, iv, c, alf, t, dummy )
c       -----------------------------------------------------
c
         do 210 j=1,n
210         bm(i,j)=ph(j)
220   do 240 j=1,npkg
         do 240 i=1,m
            sum=0.
            do 230 k=1,npkg
230            sum=sum+bm(i,k)*a(k,j)
240         xx1(i,j)=sum
      do 260 i=1,m
         sum=0.
         do 250 k=1,npkg
250         sum=sum+xx1(i,k)*bm(i,k)
         syh(i)=sqrt(abs(sum))
260      sr(i)=r(i)/sqrt(abs(rms/weight(i)-sum))
c
c-- Calculate expected standardized residuals ( for normal probability plot)
c
      if (m.gt.10) go to 410
      do 400 i=1,m
         pp=(3.*i-1.)/(3.*m+1.)
400      syh(i)=4.91*(pp**(.14)-(1.-pp)**(.14))
      go to 430
410   pp=2./(3.*m+1.)
      syh(1)=4.91*(pp**(.14)-(1.-pp)**(.14))
      pi8=pi/8
      pi4=pi/4
      do 420 i=2,m
         pp=(i-pi8)/(m+1.-pi4)
         syh(i)=4.91*(pp**(.14)-(1.-pp)**(.14))
420   continue
c    ------------------------
430   call ORPLOT (sr,m,ior)
c    ------------------------
      do 440 i=1,m
440      srv(i)=sr(ior(i))
c
c--  VG plot of (1) std. residuals vs. ind. var. 
c               (2) std. residuals vs. dep. var.
c               (3) obs. std. resids. vs. expected std. resids. 
c
      do 360 i=1,iv 
c
         do 355 j=1,m
            temp(j) = 0.0
  355    continue
c
c         if (ts.and.itp) call SCPLT2 ( t(1,i), sr, m, t(1,i), temp, m,
c     *                             0, 1, 'std. resids. vs. indep. var.',
c     *                                  tlabel(i), 'std. resids.' )
c
c         if (.not.ts) call SCPLT2 ( t(1,i), sr, m, t(1,i), temp, m, 1, 
c     *                               1, 'std. resids. vs. indep. var.',
c     *                                  tlabel(i), 'std. resids.' )
c
         if (ts.and.itp) call THRPLT ( t(1,i), sr, m, t(1,i), temp, m,
     *                             0, 1, 'Std Resids vs Ind Var',
     *                                 tlabel(i), 'Std. Resids.',
     *                                 eta, sr, m, eta, temp, m, 1, 1,
     *                                'Std Resids vs Dep Var',
     *                                 ylabel, 'Std. Resids.',
     *                     syh, srv, m, 1, 0, 'Normal Prob. Plot',
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )
c
         if (.not.ts) call THRPLT ( t(1,i), sr, m, t(1,i), temp, m, 1, 
     *                               1, 'Std Resids vs Ind Var',
     *                                  tlabel(i), 'Std. Resids.',
     *                               eta, sr, m, eta, temp, m, 1, 1,
     *                               'Std Resids vs Dep Var',
     *                               ylabel, 'Std. Resids.',
     *                    syh, srv, m, 1, 0, 'Normal Prob. Plot',
     *              'Expected Std. Resids.', 'Obs. Std. Resids.', ts )
c
c      call SCPLT2 ( eta, sr, m, eta, temp, m, 1, 1,
c     *             'Std. Resids. vs. Dependent Var.',
c     *             ylabel, 'Std. Resids.' )
c
c
c      call SCPLT1 ( syh, srv, m, 1, 0, 'Normal Probability Plot',
c     *              'Expected Std. Resids.', 'Obs. Std. Resids.' )
c
c
         if (.not.(idp.or.idcp)) go to 360
c
c--  Plot std resids. vs. ind. variable using DISSPLA graph
c
         do 350 j=1,m
            temp(j)=t(j,i)
350      continue
c
         if(ts.and.idp)call PLOT1('Std. Resids. vs. Ind. Var.$',
     2                            tlabld(i),'Standardized Residuals$',
     3                            temp, sr,m,10.,6.,0.,0.,1)
c
         if (ts.and.idcp) call PLOT1 ('Std. Resids. vs. Ind. Var.$', 
     2                            tlabld(i),'Standardized Residuals$',
     3                                temp,sr,m,8.5,5.,6.35,1.5,1)
c
         if (.not.ts.and.idp) 
     1      call PLOT1 ('Std. Resids. vs. Ind. Variable$',
     2                  tlabld(i),
     3                  'Standardized Residuals$',temp,sr,m,8.,6.,0., 
     4                  0.,1)
c 
         if(.not.ts.and.idcp) call PLOT1('Std. Resids. vs. Ind. Var.$',       
     2                                   tlabld(i),
     3                                   'Standardized Residuals$',temp,
     4                                   sr,m,7.,5.,6.35,2.,1)
c
360   continue
c
c  Plot std. resids. vs. predicted values using DISSPLA graph
c
      if (idp) call PLOT1 ('Std. Resids. vs. Pred. Values$',
     2                     ylabld,'Standardized Residuals$',eta,sr,m,8.,
     3                     6.,0.,0.,2)
c
      if (idcp) call PLOT1 ('Std. Resids. vs. Pred. Values$',ylabld,
     2                      'Standardized Residuals$',eta,sr,m,7.,5., 
     3                      6.35,2.,2)
c
c--  Do normal probability plot using DISSPLA graph
c
      if (idp) call PLOT1 ('Normal Probability Plot$',
     2                     'Expected Standardized Residuals$',
     3                     'Standardized Residuals$',syh,srv,m,8.,6., 
     4                     0.,0.,3)
c
      if (idcp) call PLOT1 ('Normal Probability Plot$',
     2                      'Expected Std. Resids.$',
     3                      'Observed Std. Resids.$',syh,srv,m,7.,5.,
     4                      6.35,2.,3)
c
c-------------------------------------------------------------------------------
c
c--  Print table of obs. & pred. values along with std. res. if desired
c
      print 40
40    format (' Print table of obs, pred val, resid & std. resid?  ')
      read (5,20) ians
20    format ( a1 )
      if ( (ians .ne. 'y') .and. (ians .ne. 'Y') ) go to 330 
c
      do 300 inv=1,iv
         print 270, inv
         if (ip) write (10,270) inv
         do 290 i=1,m
            print 280, i,t(i,inv),y(i),eta(i),r(i),sr(i)
            if (ip) write (10,280) i,t(i,inv),y(i),eta(i),r(i),sr(i)
290      continue 
300   continue
270   format ('1'/' index',2x,'ind var', i2, 4x,'obs val',5x,'pred val',
     2             6x, 'resid',8x,'std res'/)
280   format (' ',i3,1p5e13.5)
c
c-------------------------------------------------------------------------------
 330  continue
c
460   return
c
      end 
