I have a problem in estimating parameters of models. I use program r to solve these but the code that I wrote does not work. Could you help me?
This is the code that I wrote:
getwd()
setwd("D:\document\KKU Ph.D\R")
getwd()
mydata<-read.csv("dat_claims4var.csv",header = TRUE, sep=",")
attach(mydata)
#convert sex to 0,1
client_sex_n<-rep(NA,length(client_sex))
for(i in 1:length(client_sex)){
if(client_sex[i]=="Woman"){
client_sex_n[i]<-1
}else{
client_sex_n[i]<-0
}
}
client_sex_n[1]*2
#convert character to numeric
cities2_n<-as.numeric(cities2)
cities2_n[1]*2
north_n<-as.numeric(north)
north_n[1]*2
client_age_n<-as.numeric(client_age)
client_age_n[1]*2
client_nother_n<-as.numeric(client_nother)
client_nother_n[1]*2
exposi_auto_n<-as.numeric(exposi_auto)
exposi_auto_n[1]*2
nclaims_bi_n<-as.numeric(nclaims_bi)
nclaims_bi_n[1]*2
a_1<-0.5 #initial value of alpha
xbeta<-c(1,client_sex_n[1],cities2_n[1],north_n[1],client_age_n[1],client_nother_n[1])%*%(beta)
xbeta
#estimate beta_j (marginal distribution), j=1,2,3
m1.ll<-function(theta){
beta_01<-theta[1]
beta_11<-theta[2]
beta_21<-theta[3]
beta_31<-theta[4]
beta_41<-theta[5]
beta_51<-theta[6]
b<-c(beta_01,beta_11,beta_21,beta_31,beta_41,beta_51)
lxb<-c()
ll<-c()
l<-c()
for(i in 1:length(nclaims_bi_n)){
lxb[i]<-log(exposi_auto_n[i])+c(1,client_sex_n[i],cities2_n[i],north_n[i]
,client_age_n[i],client_nother_n[i])%*%b
l[i]<-log((gamma(a_1+nclaims_bi_n[i])a_1^a_1exp(nclaims_bi_n[i]*lxb[i])
/(factorial(nclaims_bi_n[i])gamma(a_1)(a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i]))))
}
ll<-sum(l,na.rm=TRUE)
return(-ll)
}
#Gradient of m1.ll
gr<-function(theta){
beta_01<-theta[1]
beta_11<-theta[2]
beta_21<-theta[3]
beta_31<-theta[4]
beta_41<-theta[5]
beta_51<-theta[6]
lxb<-c()
b<-c(beta_01,beta_11,beta_21,beta_31,beta_41,beta_51)
for(i in 1:length(nclaims_bi_n)){
#lxbb=ln(E_ij+Xbeta)
lxb[i]<-log(exposi_auto_n[i])+c(1,client_sex_n[i],cities2_n[i],north_n[i]
,client_age_n[i],client_nother_n[i])%*%b
g1<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])lxb[i])(a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
g2<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i]*client_sex_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])*lxb[i])client_sex_n[i](a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
g3<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i]*cities2_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])*lxb[i])cities2_n[i](a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
g4<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i]*north_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])*lxb[i])north_n[i](a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
g5<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i]*client_age_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])*lxb[i])client_age_n[i](a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
g6<-((factorial(nclaims_bi_n[i])*gamma(a_1)gamma(a_1+nclaims_bi_n[i])(a_1^a_1)*exp(nclaims_bi_n[i]lxb[i]))(((a_1+exp(lxb[i])^(a_1+nclaims_bi_n[i]))*nclaims_bi_n[i]*client_nother_n[i])
-(a_1+exp((a_1+nclaims_bi_n[i])*lxb[i])client_nother_n[i](a_1+exp(lxb[i])))))/(factorial(nclaims_bi_n[i])gamma(a_1)((a_1+exp(lxb[i]))^(a_1+nclaims_bi_n[i])))^2
}
c(g1, g2, g3, g4, g5, g6)
}
p0<-c(2.9,-0.1,0.1,-0.1,-0.02,-0.07)
t<-m1.ll(p0)
t
m<-optim(par=p0, m1.ll, gr, method ="L-BFGS-B" )