payments <- matrix(c(2000, 1800, 1500, 1600, 2200, 1600, 1400, 1400, 2500, 2000, 1700, 1600), nrow = 3, ncol = 4,byrow=TRUE)
payments 
eps <- 0.01

chiSquare <- function(riskFactors) {
	cartype <- c(1,riskFactors[1:2]) #c(chi11,chi12,chi13)
	age <- riskFactors[3:6] #c(chi21,chi22,chi23,chi24)
	chiSq <- 0
	for (i in 1:3){
	for (j in 1:4){
	chiSq <- chiSq + (payments[i,j] - cartype[i]*age[j])^2/(cartype[i]*age[j])
	}
	}
	return(chiSq )
}

chiSquare(c(2,2,2,2,2,2))

OptRiskFactors <- optim(c(1,1,1,1,1,1),chiSquare, lower = c(eps ,eps ,eps ,eps ,eps ,eps))$par
OptRiskFactors <- c(1,OptRiskFactors)
OptRiskFactors

tariffs <- matrix(0, nrow = 3, ncol = 4)
for (i in 1:3){
for (j in 1:4){
	tariffs[i,j] <- OptRiskFactors[i]*OptRiskFactors[3+j]
	}
	}
tariffs
payments 

totalPayments<-0
totalTariffs<-0
for (i in 1:3){
for (j in 1:4){
	totalPayments <- totalPayments + payments[i,j]
	totalTariffs <- totalTariffs + tariffs[i,j]
	}
	}
totalPayments
totalTariffs  


### Load the observed claim amounts into a matrix
 S <- matrix(c(2000 ,2200 ,2500 ,1800 ,1600 ,2000 ,1500 ,1400 ,1700 ,1600 ,1400 ,1600) , nrow =3)

 ### Define the design matrix Z
 Z <- matrix(c(rep(1,12),rep(0,4),rep(1,4),rep(0,12),rep(1,4),rep(c(0,1,0,0),3),
 rep(c(0,0,1,0),3),rep(c(0,0,0,1),3)), nrow =12)

 ### Store design matrix Z and log(S_{i,j}) in one dataset
 data <- as.data.frame(cbind(Z[,-1],matrix(log(t(S)),nrow =12)))
 colnames(data) <- c("van", "truck", "X31_40y", "X41_50y", "X51_60y", "observation")

 ### Apply the regression model
 linear.model1 <- lm(formula = observation ~ van+truck+X31_40y+X41_50y+X51_60y , data=data)
 summary(linear.model1)

 ### Fitted values
 matrix(exp(fitted(linear.model1)), byrow=TRUE , nrow =3)

 ### We can also get the parameters by applying the closed formula
 solve(t(Z)%*%Z) %*% t(Z) %*% matrix(log(t(S)), nrow =12)

 ### We can also use R directly on the data (it finds the design matrix internally)
 car <- c(" passenger car", "van", "truck")
 age <- c(" X21_30y", "X31_40y", "X41_50y", "X51_60y ")
 dat <- expand.grid(car , age)
 colnames(dat) <- c("car","age")
 dat$observation <- as.vector(log(S))
 linear.model1.direct <- lm(formula = observation ~ car+age , data=dat)
 summary(linear.model1.direct)