Difference equation approximation OU simulation function

For reasons I won't get into in detail, I just wrote a simple function to conduct Ornstein-Uhlenbeck simulations by finely discretizing time and then using a simple difference equation version of the OU model to simulate as follows:

multiOU<-function(tree,alpha,sig2,theta=NULL,a0=NULL,internal=TRUE,...){
    if(hasArg(dt)) dt<-list(...)$dt
    else dt<-1/1000*max(nodeHeights(tree))
    ss<-sort(unique(c(getStates(tree,"tips"),getStates(tree,"nodes"))))
    if(is.null(theta)) theta<-setNames(rep(0,length(ss)),ss)
    if(is.null(a0)) a0<-0
    tree<-reorder(tree,"cladewise")
    root<-Ntip(tree)+1
    S<-matrix(NA,nrow(tree$edge),2)
    S[which(tree$edge[,1]==root),1]<-a0
    for(i in 1:nrow(tree$edge)){
        x1<-S[i,1]
        for(j in 1:length(tree$maps[[i]])){
            t<-tree$maps[[i]][j]
            ALPHA<-alpha[names(t)]
            SIG2<-sig2[names(t)]
            THETA<-theta[names(t)]
            t<-c(rep(dt,floor(t/dt)),t%%dt)
            for(k in 1:length(t))
                x1<-x1+ALPHA*(THETA-x1)*t[k]+
                    rnorm(n=1,sd=sqrt(SIG2*t[k]))
        }
        S[i,2]<-x1
        if(any(tree$edge[,1]==tree$edge[i,2]))
            S[which(tree$edge[,1]==tree$edge[i,2]),1]<-S[i,2]
    }
    x<-setNames(c(S[1,1],S[,2]),c(tree$edge[1,1],
        tree$edge[,2]))[as.character(1:(Ntip(tree)+tree$Nnode))]
    names(x)[1:Ntip(tree)]<-tree$tip.label
    if(!internal) x<-x[tree$tip.label]
    x
}

In the following, I'll use it to simulate:

library(phytools)
set.seed(116)
cols<-setNames(c("black","blue","red"),letters[1:3])
tree<-pbtree(n=100,scale=2)
Q<-matrix(c(-1,1,0,1,-2,1,0,1,-1),3,3,
    dimnames=list(letters[1:3],letters[1:3]))
tree<-sim.history(tree,Q,anc="a")

## Done simulation(s).

sig2<-setNames(c(1,2,1),letters[1:3])
alpha<-setNames(c(2,2,2),letters[1:3])
theta<-setNames(c(10,0,-10),letters[1:3])
x<-multiOU(tree,alpha,sig2,theta)
phenogram(tree,x,ftype="off",colors=cols,ylim=c(-11,11),
    xlim=c(0,2))
abline(h=c(-10,0,10),lty="dashed",col=sapply(c("red","blue",
    "black"),make.transparent,alpha=0.5),lwd=3,lend=2)
add.simmap.legend(colors=cols,prompt=FALSE,x=0,y=-5)

I'm going to add (importantly) that the same simulation is done without the discrete difference equation approximation by OUwie.sim in the powerful OUwie package by Jeremy Beaulieu & Brian O'Meara.

Finally, let's fit the appropriate 'OUwie' (that is, flexible OU) model to these data:

library(OUwie)
data<-data.frame(Genus_species=tree$tip.label,
    Reg=getStates(tree,"tips"),X=as.numeric(x[tree$tip.label]))
OUwie(tree,data,model="OUMV",simmap.tree=TRUE)

## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.

## 
## Fit
##       lnL      AIC     AICc model ntax
##  -153.169 320.3379 321.5553  OUMV  100
## 
## 
## Rates
##                  a         b         c
## alpha    1.6410055  1.641006  1.641006
## sigma.sq 0.8331012 12.695801 73.799702
## 
## Optima
##                  a          b         c
## estimate 9.5391941 -2.5029839 -14.22271
## se       0.1269122  0.8906037   2.51031
## 
## Arrived at a reliable solution

One of our σ² estimates is way off - but aside from that, we seem to have more or less captured the generating process. Cool.