Hello,
I am trying to do a fit a loglikelihood function with Multivariate
distribution via copulas with fitMdvc. The problem is that it
doesn't recognize that my beta is a vector of km parameter and when I try
to run it it say that the length of my initial values is not the same as
the parameter.
Can somebody guide me where my mistake is.
Thanks,
Elisa.
#################################
#### OLS ESTIMATION
fit1 <- lm(Y1 ~ Xi)
fit2 <- lm(Y2 ~ Xi)
fit3 <- lm(Y3~ Xi)
resid1 <- Y1-as.matrix(cbind( 1,Xi))%*%fit1$coef
resid2 <- Y2-as.matrix(cbind( 1,Xi))%*%fit2$coef
resid3 <- Y3-as.matrix(cbind( 1,Xi))%*%fit3$coef
sigma1<-sum(resid1*resid1)/nrow(Xi)
sigma2<-sum(resid2*resid2)/nrow(Xi)
sigma3<-sum(resid3*resid3)/nrow(Xi)
rho12<- cor(resid1, resid2)
rho13<- cor(resid1, resid3)
rho23<- cor(resid2, resid3)
library("copula")
library(mvtnorm)
x<-as.matrix(cbind( 1,Xi))
km<- as.numeric(ncol(x))
#Parameters
beta <- rbind(rep(0,km))
sigma <- 1
rho <- 0.5
#user-defined functions
dtobit1 <- function(beta,sigma,x,y1) {ifelse(Y1>0, dnorm(Y1,x%*%beta,
sigma),(1-pnorm((x%*%beta)/sigma)))}
ptobit1 <- function(beta,sigma,x,y1) {ifelse(Y1>0,
pnorm((Y1-x%*%beta)/sigma),(1-pnorm((x%*%beta)/sigma)))}
dtobit2 <- function(beta,sigma,x,y2) {ifelse(Y2>0, dnorm(Y2,x%*%beta,
sigma),(1-pnorm((x%*%beta)/sigma)))}
ptobit2 <- function(beta,sigma,x,y2) {ifelse(Y2>0,
pnorm((Y2-x%*%beta)/sigma),(1-pnorm((x%*%beta)/sigma)))}
dtobit3 <- function(beta,sigma,x,y3) {ifelse(Y3>0, dnorm(Y3,x%*%beta,
sigma),(1-pnorm((x%*%beta)/sigma)))}
ptobit3 <- function(beta,sigma,x,y3) {ifelse(Y3>0,
pnorm((Y3-x%*%beta)/sigma),(1-pnorm((x%*%beta)/sigma)))}
#mvdc object
myMvdc <- mvdc(copula = ellipCopula("normal", param = c(rho12,
rho13,
rho23),dim = 3, dispstr = "un"),
margins = c("tobit1", "tobit2", "tobit3"),
paramMargins = list(list(beta=fit1$coef,sigma=sigma1),
list(beta=fit2$coef,sigma=sigma2), list(beta=fit3$coef,sigma=sigma3)))
start <- as.vector(c(fit1$coef,sigma1, fit2$coef,sigma2,
fit3$coef,sigma3,rho12,rho13,rho23))
#MLE
fit <- fitMvdc(as.matrix(cbind(Y1,Y2,Y3)), myMvdc , start = start ,
optim.control = list(trace = TRUE, maxit = 2000))
length(myMvdc@paramMargins)
[[alternative HTML version deleted]]
