lorelogram <- function(qdata){

#
m  <- dim(qdata)[1]; # number of individuals
n  <- dim(qdata)[2]; # number of observations (per individual)

# Values at time differences u=1,...,n-1
lor <- rep(0,times=n-1); 

# Numbers of transitions
n00    <- rep(0,times=n-1); 
n01    <- rep(0,times=n-1); 
n10    <- rep(0,times=n-1); 
n11    <- rep(0,times=n-1); 

# Loop over all time differences
for (h in 1:(n-1)){

    for (k in (h+1):n){ 
     
      # Subsets of the data according to whether
      # the h'th sample was 0 or 1
      q0  <- qdata[which(qdata[,h]==0),];
      q1  <- qdata[which(qdata[,h]==1),];


      # Differece between the k'th and h'th samples
      apu0 <- q0[,k]-q0[,h];
      apu1 <- q1[,k]-q1[,h];

      # Update the numbers of transitions
      n00[k-h] <- n00[k-h] + sum(apu0==0);
      n01[k-h] <- n01[k-h] + sum(apu0==1);
      n10[k-h] <- n10[k-h] + sum(apu1==-1);
      n11[k-h] <- n11[k-h] + sum(apu1==0);
   }  
    
}
 
# Calculate the odds ratios for each time difference
for (i in 1:(n-1)){
    lor[i] <- (n11[i]/n10[i])/(n01[i]/n00[i]);
}

# Plot the log odds ratio vs. time difference
plot(log(lor));
lines(log(lor));

# Return the values of the log odds ratios
return(log(lor))

}