R does not display output


#1

Am trying to run this code but R does nor display any output , the code is correct but no output ?

BerryLevinsohnPakes <- function(dat, mkt.id.fld = "mkt.id", prod.id.fld = "prod.id", prc.fld = "px", share.fld = "share", x.var.flds = c("x1", "x2", "x3"), prc.iv.flds = "z", n.sim = 1500, sigma.guess){
# all data should appear in the data.frame "dat"
# input variables ending in ".fld" are the names of the columns
# n.sim = number of simulated "indviduals" per market
# sigma guess may be missing
#
# Required packages (installed as the function runs):
# SQUAREM, BB, AER
#
# Based on code written by Aviv Nevo, May 1998.
# Adapted by Michael Carniol, January 2015

blp_inner <- function(delta.in, mu.in) {
	# Computes a single update of the BLP (1995) contraction mapping.
    # of market level predicted shares.
	# This single-update function is required by SQUAREM, see Varadhan and
	# Roland (SJS, 2008), and Roland and Varadhan (ANM, 2005)
	# INPUT
	# 	delta.in : current value of delta vector
	# 	mu.in: current mu matrix
	# Requires global variables: s.jt
	# OUTPUT
	# 	delta.out : delta vector that equates observed with predicted market shares
    pred.s <- rowMeans(ind_sh(delta.in, mu.in));
    delta.out <- delta.in + log(s.jt) - log(pred.s)
    return(delta.out)
}
ind_sh <- function(delta.in, mu.in){
	# This function computes the "individual" probabilities of choosing each brand
	# Requires global variables: mkt.id, X, v
	numer <- exp(mu.in) * matrix(rep(exp(delta.in), n.sim), ncol = n.sim);
	denom <- as.matrix(do.call("rbind", lapply(mkt.id, function(tt){
		1 + colSums(numer[mkt.id %in% tt, ])
	})))
	return(numer / denom);	
}
gmm_obj <- function(theta2){
	# This function computes the GMM objective function
	# Requires global variable inputs: X, v, delta, a, W
	# Outputs: theta1, xi.hat
	print(paste0("GMM Loop number: ", Sys.time()))
	print(a <<- a + 1)
	print("Updated theta2 estimate:")
	print(theta2)
	print("Change in theta2 estimate:")
	print(theta.chg <- as.numeric(theta2 - theta2.prev));
	if(sum(theta.chg != 0) <= 2){
		delta <- dat[, "delta"];
	} else {
		delta <- Y;
	}
	theta2.prev <<- theta2;

    mu <- X %*% diag(theta2) %*% v;

	print("Running SQUAREM contraction mapping")
	print(system.time(
		squarem.output <- squarem(par = delta, fixptfn = blp_inner, mu.in = mu, control = list(trace = TRUE))
	));
	delta <- squarem.output$par
	print(summary(dat[, "delta"] - delta));
	dat[, "delta"] <<- delta;
	
	mo.ivreg <- ivreg(fm.ivreg, data = dat, x = TRUE)
	theta1 <<- coef(mo.ivreg);
	xi.hat <<- as.vector(mo.ivreg$resid);
	Z.hat <- Z * matrix(rep(xi.hat, ncol(Z)), ncol = ncol(Z))
	W.inv <- try(solve(t(Z.hat) %*% Z.hat), silent = FALSE)
	if("matrix" == class(W.inv)){
		PZ <<- Z %*% W.inv %*% t(Z);
		PX.inv <- solve(t(X) %*% PZ %*% X)
		theta1 <<- PX.inv %*% t(X) %*% PZ %*% delta
		xi.hat <<- delta - X %*% theta1
		X.hat <- (PZ %*% X) * matrix(rep(xi.hat, K), ncol = K)
		tsls.se <- sqrt(diag(PX.inv %*% t(X.hat) %*% X.hat %*% PX.inv))
		print("GMM step 2 updated theta1 estimate:")
		print(beta.est <<- data.frame(beta.est = theta1, beta.se = tsls.se, sigma.est = theta2))
	}
	dat[, "xi.hat"] <<- xi.hat
	f <- t(xi.hat) %*% PZ %*% xi.hat;
	print("Updated GMM objective:")
	print(f <- as.numeric(f));
	return(f)
}
jacobian <- function(delta.in, theta.in){
	print(paste0("Calculating Jacobian matrix, ", Sys.time()))
	#Requires global variables X, v, mkt.id
	mu1 <- X %*% diag(theta.in) %*% v;
	ind.shares <- ind_sh(delta.in, mu1);
	K <- ncol(X);
	print(paste0("Calculating dsigma matrix, ", Sys.time()))
	dsigma <- lapply(l.Xv, function(x){
		temp2 <- x * ind.shares;
		temp3 <- as.matrix(do.call("rbind", lapply(mkt.id, function(m){
			colSums(temp2[mkt.id %in% m, ])
		})));
		dsigma.res <- rowMeans(temp2 - ind.shares * temp3);
		return(dsigma.res)
	})
	dsigma <- as.matrix(do.call("cbind", dsigma))
	print(paste0("Calculating ddelta matrices, ", Sys.time()))
	ddelta <- list()
	for(m in mkt.id){
		if(m %in% names(ddelta)){next}
		temp1 <- as.matrix(ind.shares[mkt.id %in% m, ]);
		H1 <- temp1 %*% t(temp1);
		H2 <- diag(rowSums(temp1));
		H <- (H2 - H1) / n.sim;
		H.inv <- solve(H);
		ddelta[[as.character(m)]] <- H.inv %*% dsigma[mkt.id %in% m, ];
		rm(temp1, H1, H2, H, H.inv)
	}
	ddelta <- as.matrix(do.call("rbind", ddelta));
	return(ddelta)
}
gradient_obj <- function(theta2){
	#Requires global variables PZ, delta, xi.hat
	print(system.time(jacobian_res <<- jacobian(as.vector(dat[, "delta"]), theta2)))
	print(paste0("Updated gradient:", Sys.time()))
	print(f <- -2 * as.numeric(t(jacobian_res) %*% PZ %*% xi.hat));
	#######
	L <- ncol(Z)
	covg <- matrix(0, nrow = L, ncol = L)
	for(i in 1:JT){
		covg <- covg + (Z[i, ] %*% t(Z[i, ])) * xi.hat[i]^2
	}
	d.delta <- jacobian_res;
	Dg <- t(d.delta) %*% Z
	p.Dg <- try(solve(Dg %*% W.inv %*% t(Dg)))
	cov.mat <- p.Dg %*% (Dg %*% W.inv %*% covg %*% W.inv %*% t(Dg)) %*% p.Dg
	beta.est$sigma.se <<- sqrt(diag(cov.mat));
	print(paste0("Updated coefficients table:", Sys.time()))
	print(beta.est)
	write.csv(beta.est, file = paste0("BLP_beta_est_", Sys.Date(), ".csv"))
	#######
	return(as.numeric(f))
}
	
#Set up data
dat <- dat[dat[, share.fld] > 0, ]
dat <- dat[order(dat[, mkt.id.fld], dat[, prod.id.fld]), ]
JT <- nrow(dat)
#market identifier variable
mkt.id <- dat[, mkt.id.fld];
#Number of characteristics (including constant and price)
X <- as.matrix(cbind(ones = rep(1, JT), dat[, c(x.var.flds, prc.fld)]));
K <- ncol(X)
#Compute the outside good market share by market
s.jt <- as.vector(dat[, share.fld]);
temp <- aggregate(s.jt, by = list(mkt.id = mkt.id), sum);
sum1 <- temp$x[match(mkt.id, temp$mkt.id)];
s.j0 <- as.vector(1 - sum1);
rm(temp, sum1);
dat[, "delta"] <- Y <- log(s.jt) - log(s.j0);
iv <- dat[, prc.iv.flds]

while(!require(AER)){install.packages("AER")}
#Construct 2SLS regression specification
str.ivreg.y <- "delta ~ "
str.ivreg.x <- paste(x.var.flds, collapse = " + ")
str.ivreg.prc <- paste(prc.fld, collapse = " + ")
str.ivreg.iv <- paste(prc.iv.flds, collapse = " + ")
print("2SLS specification:")
print(fm.ivreg <- paste0(str.ivreg.y, str.ivreg.x, " + ", str.ivreg.prc, " | ", str.ivreg.x, " + ", str.ivreg.iv))
rm(str.ivreg.y, str.ivreg.x, str.ivreg.prc, str.ivreg.iv)

print("2SLS beta estimate:")
print(summary(mo.ivreg <- ivreg(fm.ivreg, data = dat, x = TRUE)))
beta.est <- summary(mo.ivreg)$coef[, 1:2]
#Z = instrumental variable matrix include exogenous X's
Z <- as.matrix(mo.ivreg$x$instruments)
PZ <- Z %*% solve(t(Z) %*% Z) %*% t(Z);
theta1 <- coef(mo.ivreg);
xi.hat <- as.vector(mo.ivreg$resid);
Z.hat <- Z * matrix(rep(xi.hat, ncol(Z)), ncol = ncol(Z))
W.inv <- try(solve(t(Z.hat) %*% Z.hat), silent = FALSE)
if("matrix" == class(W.inv)){
	PZ <- Z %*% W.inv %*% t(Z);
	PX.inv <- solve(t(X) %*% PZ %*% X)
	theta1 <- PX.inv %*% t(X) %*% PZ %*% Y
	xi.hat <- Y - X %*% theta1
	X.hat <- (PZ %*% X) * matrix(rep(xi.hat, K), ncol = K)
	tsls.se <- sqrt(diag(PX.inv %*% t(X.hat) %*% X.hat %*% PX.inv))
	print("GMM step 2 updated theta1 estimate:")
	print(beta.est <- data.frame(beta.est = theta1, se.est = tsls.se))
}
dat[, "xi.hat"] <- xi.hat

#Starting point
print("Sigma guess:")
if(missing(sigma.guess)){
	tsls.se <- beta.est[, 2]
	print(theta2 <- 0.5 * tsls.se);
} else {
	print(theta2 <- sigma.guess);
}
theta2.prev <- theta2;

# Matrix of individuals' characteristics
#	Standard normal distribution draws, one for each characteristic
v <- matrix(rnorm(K * n.sim), nrow = K, ncol = n.sim)
# Break X and v matrices into list variables 
# in attempt to expedite calculation of the Jacobian matrix
l.X <- lapply(1:K, function(k){
	return(X[, k])
})
l.v <- lapply(1:K, function(k){
	return(v[k, ])
})
l.Xv <- lapply(1:K, function(k){
	l.X[[k]] %*% t(l.v[[k]]);
})

print("Estimating random coefficients multinomial logit")
a <- 0;
beta.est <- NULL;
while(!require(SQUAREM)){install.packages("SQUAREM")}
while(!require(BB)){install.packages("BB")}
print(system.time(
	theta.est <- multiStart(par = theta2, fn = gmm_obj, gr = gradient_obj, lower = 0, control = list(trace = TRUE), action = "optimize")
));
save(theta.est, file = paste0("theta_est_", Sys.time(), ".RData"))

print("Final coefficients estimate:")
theta2 <- theta.est$par
gmm.res <- gmm_obj(theta2)
grad.res <- gradient_obj(theta2)	
return(list(coef.mat = beta.est, gmm.obj.func = gmm.res, gmm_est = theta.est, final.data = dat))

}


#2

Could you please turn this into a self-contained reprex (short for minimal reproducible example)? It will help us help you if we can be sure we're all working with/looking at the same stuff.

Right now the best way to install reprex is:

# install.packages("devtools")
devtools::install_github("tidyverse/reprex")

If you've never heard of a reprex before, you might want to start by reading the tidyverse.org help page. The reprex dos and don'ts are also useful.

If you run into problems with access to your clipboard, you can specify an outfile for the reprex, and then copy and paste the contents into the forum.

reprex::reprex(input = "fruits_stringdist.R", outfile = "fruits_stringdist.md")

For pointers specific to the community site, check out the reprex FAQ, linked to below.