Я хотел бы запустить симуляцию по методу Монте-Карло в 2000 раз по статистике правдоподобия c для следующей системы уравнений:
y1 = a1 + b1X + c1K
y2 = a2 + b2X + c2K
Я хочу применить ограничение. Гипотеза, которую я хочу проверить: Ho: a1 + 39 = a2.
Предположим, что каждая выборка имеет 50 наблюдений.
Я использовал пакет systemfit для расчета моей имитированной модели SUR. Я просто не уверен, как применить ограничение в R.
Вот как я начал моделировать модель с ограничениями:
#Set Parameters for Monte Carlo Simulation
N <- 50
nrsim <- 2000
#Equation 1
a1 <- -29
b1 <- 0.12
g1 <- 0.12
mu1 <- 0
sig2e1 <- 1700
sde1 <- sqrt(sig2e1)
F1 <- unlist(mydata %>%
filter(eqn == 1) %>%
select(`F`))
C1 <- unlist(mydata %>%
filter(eqn == 1) %>%
select(C))
va1 <- numeric(nrsim) #stores the estimates of alpha1
vb1 <- numeric(nrsim) #stores the estimates of beta1
vg1 <- numeric(nrsim) #stores the estimates of gamma1
#Equation 2
a2 <- -12
b2 <- 0.078
g2 <- 0.9
mu2 <- 0
sig2e2 <- 1900
sde2 <- sqrt(sig2e2)
F2 <- unlist(mydata %>%
filter(eqn == 2) %>%
select(`F`))
C2 <- unlist(mydata %>%
filter(eqn == 2) %>%
select(C))
#unrestricted model
unr.vsura1 <- numeric(nrsim) #stores the estimates of alpha1
unr.vsura2 <- numeric(nrsim) #stores the estimates of alpha2
unr.vsurb1 <- numeric(nrsim) #stores the estimates of beta1
unr.vsurb2 <- numeric(nrsim) #stores the estimates of beta2
unr.vsurg1 <- numeric(nrsim) #stores the estimates of gamma1
unr.vsurg2 <- numeric(nrsim) #stores the estimates of gamma2
for (j in 1:nrsim) {
set.seed(12345+10*j)
y1 <- as.numeric(unlist(a1 + b1*F1 + g1*C1 + rnorm(N, mean = mu1, sd = sde1)))
y2 <- as.numeric(unlist(a2 + b2*F2 + g2*C2 + rnorm(N, mean = mu2, sd = sde2)))
X1 <- cbind(F1, C1)
X2 <- cbind(F2, C2)
eq1 <- y1~X1
eq2 <- y2~X2
system <- list(eq1 = eq1, eq2 = eq2)
#model
sur <- systemfit(system, method = "SUR", data = mydata)
#store each iteration of the estimate
unr.vsura1[j] <- coef(sur)[[1]]
unr.vsura2[j] <- coef(sur)[[4]]
unr.vsurb1[j] <- coef(sur)[[2]]
unr.vsurb2[j] <- coef(sur)[[5]]
unr.vsurg1[j] <- coef(sur)[[3]]
unr.vsurg2[j] <- coef(sur)[[6]]
}