rm(list=ls())

library(ape)
library(geiger)
library(extraDistr)
library(LaplacesDemon)
#install.packages("RRphylo")
library(RRphylo)


ntaxa<-3


kappa<-0.002
alpha<-0.006
tausq<-1
betas<-rnorm(2*ntaxa-1)

tree<-rcoal(ntaxa)

plot(tree)
nodelabels()
tiplabels()
edgelabels(round(tree$edge.length,2))


makeL(tree)->L
round(L,2)


makeL1(tree)->L1
round(L1,2)


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]) 
  # + alpha/treelength[index]*betas[anc[index]]
}  

trait


# up the function for parameter estimation for ar1, then do testing

#Hessian for like ar(1) testing

# Ar1 can be test two rate wood and plant tree seems cool