c ****************************************************************** c OLS Regression c ****************************************************************** subroutine olsx(y,x,n,betahat,rmse,designx,xpx,xpy,xpxinv) implicit real*8(a-h,o-z) real*8 y(n,1),x(n,1),betahat(2,1),designx(n,2) real*8 xpx(2,2),xpy(2,1),xpxinv(2,2) call zeros(designx,n,2) call zeros(xpx,2,2) call zeros(xpy,2,1) do 90 i=1,n designx(i,1) = 1 designx(i,2) = x(i,1) 90 continue do 200 i=1,n do 150 j=1,2 xpy(j,1) = xpy(j,1) + (designx(i,j) * y(i,1)) 150 continue 200 continue do 300 i=1,n do 250 j=1,2 do 225 k=1,2 xpx(j,k) = xpx(j,k) + (designx(i,j) * designx(i,k)) 225 continue 250 continue 300 continue call xinverse(xpx,xpxinv,2,bigdet) call zmmultm(xpxinv,xpy,betahat,2,2,1) rmse = 0 do 400 i=1,n rresid = y(i,1) - betahat(1,1) - (betahat(2,1) * x(i,1)) rmse = rmse + (rresid ** 2) 400 continue rmse = dsqrt(rmse / (n - 2)) return end c ****************************************************************** c Method of moments estimate of sigma_u c ****************************************************************** double precision function findsigu(ww,msamp,numw) implicit real*8(a-h,o-z) real*8 ww(msamp,numw) sigsqu = 0 do 100 i=1,msamp zmu = 0 do 10 j=1,numw zmu = zmu + (ww(i,j) / dfloat(numw)) 10 continue do 20 j=1,numw sigsqu = sigsqu + (( (ww(i,j) - zmu) ** 2) 6 / (msamp * (numw - 1))) 20 continue 100 continue findsigu = dsqrt(sigsqu) return end c ***************************************************************** c Method of moments estimate accounting for sampling weights c ***************************************************************** subroutine mom(yy,ww,siguest,pweight,n,msamp,numw,dummyw, 6 sigeest,beta0,beta1,sigxmom,zmuxmom) implicit real*8(a-h,o-z) real*8 yy(n,1),ww(msamp,numw),pweight(n,1), 6 dummyw(msamp,numw) c Compute the sum of the weights for those sampled pwsum = 0 do 10 i=1,msamp 10 pwsum = pwsum + pweight(i,1) c Compute the mean of Y. Does not need weights. zmuy = zmeanc(yy,n,1,1) print *, "Mean of Y = ",zmuy c Compute the s.d. of Y. Does not need weights. zsigy = stdc(yy,n,1,1) print *, "s.d. of Y = ",zsigy c Compute the mean of W. This does need weights zmuw = 0 do 100 i=1,msamp dummyw(i,1) = 0 do 90 j=1,numw zmuw = zmuw + (pweight(i,1) * ww(i,j) 6 / (pwsum * numw)) 90 dummyw(i,1) = dummyw(i,1) + (ww(i,j) / numw) 100 continue zmuxmom = zmuw print *, "Moments estimate of mean of W = ",zmuw c Compute the s.d. of Wbar. This must account for the weights. zsigwbar = 0. do 150 i=1,msamp 150 zsigwbar = zsigwbar 6 + (pweight(i,1) * ((dummyw(i,1) - zmuw) ** 2) 6 / pwsum) print *, "Moments estimate of var of W_bar = ",zsigwbar zsigwbar = dsqrt(zsigwbar) c Compute the s.d. of X zsigx = dsqrt((zsigwbar ** 2) - ( (siguest ** 2) / numw)) sigxmom = zsigx c Compute the covariance of Y and X. This must account c for the weights zsigyx = 0 do 200 i=1,msamp 200 zsigyx = zsigyx + 6 ( pweight(i,1) * 6 (yy(i,1) - zmuy) * (dummyw(i,1) - zmuw) / pwsum) beta1 = zsigyx / (zsigx ** 2) beta0 = zmuy - (beta1 * zmuw) sigesq = (zsigy ** 2) - ( (beta1 * zsigx) ** 2) sigeest = dsqrt(sigesq) return end c ***************************************************************** c Method of moments estimate accounting for sampling weights c ***************************************************************** subroutine mom2(yy,ww,siguest,pweight,n,msamp,numw,dummyw, 6 sigeest,beta0,beta1,sigxmom,zmuxmom,ierrcde) implicit real*8(a-h,o-z) real*8 yy(n,1),ww(msamp,numw),pweight(n,1), 6 dummyw(msamp,numw) ierrcde = 0 c Compute the sum of the weights for those sampled pwsum = 0. do 10 i=1,msamp 10 pwsum = pwsum + pweight(i,1) do 30 i=1,msamp dummyw(i,1) = 0. do 20 j=1,numw 20 dummyw(i,1) = dummyw(i,1) + (ww(i,j) / numw) 30 continue zmuhatx = 0. do 40 i=1,msamp 40 zmuhatx = zmuhatx + (pweight(i,1) * dummyw(i,1) / pwsum) sighxsq = 0. do 50 i=1,msamp 50 sighxsq = sighxsq + 6 ( pweight(i,1) * ((dummyw(i,1) - zmuhatx) ** 2)) sighxsq = (sighxsq / pwsum) - 6 (siguest * siguest / numw) if (sighxsq.ge.0.0d0) sigxest = dsqrt(sighxsq) if (sighxsq.lt.0.0d0) then sigxest = .01 ierrcde = 1 endif c print *, "Estimated Variance of X = ",sighxsq c print *, "Estimated s.d. of X = ",sigxest zmuxmom = zmuhatx sigxmom = sigxest a11 = 0 a12 = 0 a21 = 0 a22 = 0 b1 = 0 b2 = 0 do 60 i=1,msamp a11 = a11 + pweight(i,1) a12 = a12 + (pweight(i,1) * dummyw(i,1)) a21 = a21 + (pweight(i,1) * dummyw(i,1)) a22 = a22 + (pweight(i,1) * 6 ((dummyw(i,1) ** 2) - ((siguest ** 2)/numw))) b1 = b1 + (pweight(i,1) * yy(i,1)) 60 b2 = b2 + (pweight(i,1) * yy(i,1) * dummyw(i,1)) cdet = (a11 * a22) - (a12 * a21) c11 = a22 / cdet c22 = a11 / cdet c12 = -a12 / cdet c21 = -a21 / cdet beta0 = (c11 * b1) + (c12 * b2) beta1 = (c21 * b1) + (c22 * b2) sighsqe = 0. do 90 i=1,msamp aa = (yy(i,1) - beta0 - (beta1 * dummyw(i,1))) ** 2 90 sighsqe = sighsqe + (pweight(i,1) * aa / pwsum) sighsqe = sighsqe - (beta1 * beta1 * (siguest ** 2) / numw) if (sighsqe.ge.0.0d0) sigeest = dsqrt(sighsqe) if (sighsqe.lt.0.0d0) then ierrcde = 1 sigeest = 0.001 endif c c Small sample modifications a la Fuller c mu_y = b1 / a11 c mu_w = a12 / a11 c siguu = sigma_u^2 / numw c zmuy = b1 / a11 zmuw = a12 / a11 siguu = (siguest ** 2) / numw zmyy = 0 zmww = 0 zmyw = 0 zzz = 1.0d0 + (1.0d0 / (dfloat(msamp) - 1.0d0)) zzz2 = (1.0d0 / (dfloat(msamp) - 1.0d0)) do 110 i=1,msamp zmyy = zmyy + (pweight(i,1) * ((yy(i,1) - zmuy)**2) / a11) zmww = zmww + (pweight(i,1) * 6 ((dummyw(i,1) - zmuw)**2) / a11) zmyw = zmyw + (pweight(i,1) * (yy(i,1) - zmuy) * 6 (dummyw(i,1) - zmuw) / a11) 110 continue zlambda = ( (zmyy * zmww) - (zmyw ** 2)) / (zmyy * siguu) hhat = zmww - siguu beta1 = zmyw / hhat if (zlambda.lt.zzz) hhat=zmww - (siguu * (zlambda - zzz2)) beta1 = zmyw / (hhat + (2 * siguu * zzz2)) if (zlambda.lt.1.0d0) sigeest = 0.001 if (zlambda.lt.1.0d0) sigxmom = zmww - (zlambda * siguu) return end