rm(list=ls())

library(ape)
library(geiger)
library(extraDistr)
library(LaplacesDemon)

ntaxa<-5
kappa<-0.002
alpha<-0.006
tausq<-1
betas<-rnorm(2*ntaxa-1)
tree<-rtree(ntaxa)
makeL(tree)->L
makeL1(tree)->L1
names(betas)<-colnames(L)
trait<-array(NA,c(2*ntaxa-1)) #this is Garch(1,1)
names(trait)<-colnames(L)
betas<-betas[c(ntaxa:(2*ntaxa-1),1:(ntaxa-1))]
trait<-trait[c(ntaxa:(2*ntaxa-1),1:(ntaxa-1))]
trait[names(trait)==ntaxa+1] <- betas[ntaxa+1]
betas
trait

tree<-reorder(tree,"postorder")
tree$root.edge<-0
N<-dim(tree$edge)[1]
anc<-tree$edge[,1]
des<-tree$edge[,2]
treelength<-tree$edge.length


for(index in N:1){
  trait[des[index]]<- betas[des[index]]+ kappa/treelength[index]*betas[anc[index]]*rnorm(1,0,treelength[index])^2  + alpha/treelength[index]*betas[anc[index]]
}  

trait