# 
## G-estimation for general K-interval problem with linear blip
## 	Uses my recursive method of calculating the variance
## 	         (Completely randomized treatment) 
## 		Written by Erica E M Moodie, 16 June, 2005 		 
#  

g_estCov_1int_SeqRand.fcn <- function(Y,status,action,no.stat.vbls=1,plot.H=F,return.U=F,true.psi) {
   n <- length(Y)
   if(plot.H) par(mfrow=c(2,round(K/2)))
   
   n.stat <- no.stat.vbls
   L <- as.matrix(status)[,1:n.stat]
   if(!is.matrix(L)) L <- matrix(L,ncol=1)
   A <- action
   Y.dash <- Y

   p.coef <- coef(glm(A~L,family=binomial)) 
   p.A.given.L <- expit(cbind(1,L)%*%p.coef)  
   p.temp <- exp(cbind(1,L)%*%p.coef)  
   aa <- as.vector(A - p.A.given.L)

   B <- mean(Y.dash*aa)
   B <- rbind(B,mean(Y.dash*L*aa))
   X <- matrix(nrow=2,ncol=2)   
   X[1,1] <- mean(A*aa)
   X[2,1] <- X[1,2] <- mean(A*L*aa)
   X[2,2] <- mean(A*L^2*aa)

   W <- cbind(1,L)
   W2.inv <- solve(matrix(c(1,mean(L),mean(L),mean(L^2)),ncol=2))

   C1.temp <- apply(matrix(rep(cbind(mean(Y.dash),mean(Y.dash*L)) %*% solve(matrix(c(1,mean(L),mean(L),mean(L^2)),ncol=2)),each=n),ncol=2) * cbind(rep(1,n),L),1,sum)
   C1 <- rbind(mean(C1.temp*aa),mean(C1.temp*L*aa))

   C2.temp <- W2.inv %*% matrix(c(mean(A),mean(A*L),mean(A*L),mean(A*L^2)),ncol=2) 
   C2 <- matrix(rep(0,4),ncol=2)
   if(n<500) {
      for(i in 1:n) C2 <- C2 + rbind(aa[i],(L*aa)[i]) %*% cbind(1,L[i]) %*% C2.temp
      C2 <- C2/n
   } else {
      Mat1<-rbind(aa,L[,1]*aa)
      Mat2<-cbind(0*L[,1]+1,L[,1])
      FMat<-matrix(0,2,2)
      nblocks<-n/500
      for(i in 1:nblocks){
         i1<-(i-1)*500+1
         i2<-i*500
         sub1<-Mat1[,i1:i2]
         sub2<-Mat2[i1:i2,]
         FMat<-FMat+sub1 %*% (sub2 %*% C2.temp)
      }
      C2<-FMat/n
   }

   psi.hat <- solve(X-C2)%*%(B-C1)

   S.part <- cbind(aa,aa*L)
   d.alpha.S <- diag(c(mean(-p.A.given.L/(1+p.temp)),mean(-p.A.given.L/(1+p.temp)*L^2)))
   d.alpha.S[1,2] <- d.alpha.S[2,1] <- mean(-p.A.given.L/(1+p.temp)*L)

   Hj <- Y.dash - (A*W %*% psi.temp)
   Hj.ctr <- Hj - W %*% (W2.inv %*% rbind(mean(Hj),mean(Hj*L))) 

   d.alpha.U <- diag(-c(mean(Hj.ctr*p.A.given.L/(1+p.temp)),mean(Hj.ctr*p.A.given.L/(1+p.temp)*L^2)))
   d.alpha.U[1,2] <- d.alpha.U[2,1] <- -mean(Hj.ctr*p.A.given.L/(1+p.temp)*L)
   d.varsigma.U <- -t(W)%*%(W*aa)/n

   U.adj <- t(cbind(Hj.ctr*aa,Hj.ctr*aa*L)) - t(d.alpha.U/d.alpha.S) %*% t(as.vector(S.part)) + 
       d.varsigma.U %*% W2.inv %*% t(cbind(Hj.ctr,Hj.ctr*L)) 

   d.psi.U <- -t(W)%*%(W*aa*A)/n - 
       d.varsigma.U %*% W2.inv %*% t(W)%*%(W*A)/n  # Chain rule
   temp.a <- as.vector(solve(d.alpha.S) %*% apply(S.part,2,mean))
   x.term.for.alpha <- diag(c(mean(temp.a[1]*A + temp.a[2]*A*L),mean(temp.a[1]*A*L^2 + temp.a[2]*A*L^3)))
   x.term.for.alpha[1,2] <- x.term.for.alpha[2,1] <- mean(temp.a[1]*A*L + temp.a[2]*A*L^2)
   d.psi.U.adj <- d.psi.U - x.term.for.alpha + d.varsigma.U %*% W2.inv %*% t(W)%*%(W*A)/n

   var.psi <- diag((solve(d.psi.U.adj) %*% U.adj %*% t(U.adj) %*% t(solve(d.psi.U.adj)))/(n^2))

   if(plot.H & n.stat==1) plot(L,Y.dash-A*cbind(1,L)%*% psi.temp[1:2],ylab="H_j(psi)",xlab="L_j",main="j = 1")

   ## Testing U with TRUE psi's  
   if(return.U) {
      m <- length(true.psi)
      Hj.true <- Y.dash - A*(cbind(1,L) %*% true.psi[1:2])
      Hj.ctr <- Hj.true - W %*% (W2.inv %*% apply(as.vector(Hj) * W,2,mean))
      U.dagger <- t(cbind(Hj.ctr*aa,Hj.ctr*aa*L)) 
      U.hat <- apply(U.dagger,1,mean)

      d.a.u <- diag(-c(mean(Hj.ctr*p.A.given.L/(1+p.temp)),mean(Hj.ctr*p.A.given.L/(1+p.temp)*L^2)))
      d.a.u[1,2] <- d.a.u[2,1] <- -mean(Hj.ctr*p.A.given.L/(1+p.temp)*L)
      uu <- U.dagger - t(d.a.u) %*% solve(d.alpha.S) %*% t(S.part) + 
        d.varsigma.U %*% W2.inv %*% t(cbind(Hj.ctr,Hj.ctr*L)) 
      var.U.hat <- diag(uu %*% t(uu)/(n^2))

      return(list(psi.hat=psi.hat,psi.var=var.psi,U.hat=U.hat,U.var.hat=var.U.hat))
   }
   return(list(psi.hat=psi.hat,psi.var=var.psi))

}