library("statnet")

data("faux.mesa.high")
myformula <- faux.mesa.high ~ edges + nodematch("Grade", diff=TRUE)
mymodel <- ergm(myformula)
mymodel$coef

mynetwork1 <- simulate(mymodel, constraints = ~ edges)

mynetwork2 <- simulate(mymodel, constraints = ~ edges + bd(maxout = 5))

mynetwork3 <- faux.mesa.high
mynetwork3[ , ] <- 0 
for(i in 1:203) mynetwork3[i,i+2] <- mynetwork3[i+2, i] <- 1 

mynetwork4 <- simulate(mynetwork3 ~ edges + nodematch("Grade", diff=TRUE),
theta0 = mymodel$coef, constraints = ~edges + bd(maxout=5), seed = 12345)

summary(mynetwork4 ~ edges + degree(0:12))

mynetwork5 <- simulate(mynetwork3 ~ edges + nodematch("Grade", diff=TRUE),
theta0 = mymodel$coef, constraints = ~ bd(maxout=5), seed = 12345,
control = control.simulate(prop.weights = "random"), burnin = 1e+5)
summary(mynetwork5 ~ edges + degree(0:12))