i have monte carlo problem need compare normal curve, have this:
monte_carlo <- function(t){ k <- 2 beta1 <- 0.131 beta2 <- 0.406 betaa <- c(beta1,beta2) beta<-as.vector(betaa) phi <- 0.6 c <- 2 ni<- as.vector(rnorm(t)) x0 <- rnorm(1,(c/(1-phi)),(1/(1-(phi)^2)^(0.5))) # gerar x x_3 <- rep(0,t) x_3[1]<- c+phi*x0+ni[1] for(i in 2:t){ x_3[i]<-c+phi*x_3[i-1]+ni[i] } # fazer o loop n <- 100000 vest <- matrix(,100000,k) vt <- vest for(j in 1:n){ # gerando e e <- as.vector(rnorm(t)) # gerando y y<- rep(0,t) for(i in 1:t){ y[i]<- beta1+beta2*x_3[i]+e[i] } x <- matrix(1,t,k) x[,2]<- x_3 # estimando b por mqo b.est<- function(y,x){ xy <- t(x)%*%y # x'y xxi <- solve(t(x)%*%x) #(x'x)^(-1) b <- as.vector(xxi%*%xy) return(b) } # chamando os resultados vest[j,] <- b.est(y,x) } # obtendo médias b_mean <- colmeans(vest) # plotando os gráficos g_b1<- plot(density(vest[,1])) g_b2<-plot(density(vest[,2])) table <- cbind(medias_dos_estimadores=b_mean, grafico_b1=g_b1, grafico_b2=g_b2) return(list(b=table)) } and get, example, this:
> monte_carlo(23) $b medias_dos_estimadores media_das_estatisticas_t [1,] 0.1230776 -0.008260573 [2,] 0.4076724 0.001772357 
however, need plot normal curve compare. how should plot it?
ps.: sd in case given by:
var.est<- function(y,x){ x <- as.matrix(cbind(int=1,x)) y <- as.vector(y) t <- length(y) k <- ncol(x)-1 xy <- t(x)%*%y # x'y xxi <- solve(t(x)%*%x) #(x'x)^(-1) b <- as.vector(xxi%*%xy) #estimated coefficients yhat <- as.vector(x%*%b) #predicted values y res <- y-yhat #model residual sse <- t(res)%*%res # or sum(res^2) t(res)%*%res sigma2 <- sse/(t-k) s2 <- sigma2[1,1] est.var <- s2*solve(crossprod(x)) #coefficient standard errors var.est <- rep(0,k) for(i in c(1:k)){ var.est[i] <- est.var[i,i] } return(var.est) }
