Simulating data generating process correctly
The first thing that jumps out from the definition is that, if you have X, A, B and epsilon, you can compute y deterministically. This means you can readily test your implementation. You should always strive to find ways to define pure functions in your R code, and try to use vectorisation instead of for loops.

Based on your existing code, I'll assume for a given model that X is a matrix (n_sample, 100), A and B are vectors of length 100 and epsilon is a vector of length n_sample.

Based on your implementation, the function would look something like

compute_y &lt;- function(X, A, B, epsilon) {
  n_sample &lt;- nrow(X)
  # note that your f_[1|2] stored `double`s not `integers`
  f_1 &lt;- numeric(n_sample)
  f_2 &lt;- numeric(n_sample)

  for (d in seq(100)){
    part1 &lt;- 1 + (-1)^d*A[d]*X[,d] + B[d]*sin(6*X[,d])
    f_1 &lt;- f_1 + part1
  }
  for (d in seq(50)){
    part2 &lt;- 1 + X[,d]/50
    f_2 &lt;- f_2 + part2
  }

  f_1 * f_2 + epsilon
}


But that's a bit scruffy.

The easiest bit to clean up is the bit that defines f_2:

f_2 &lt;- numeric(n_sample)
for (d in seq(50)) {
    part2 &lt;- 1 + X[,d]/50
    f_2 &lt;- f_2 + part2
}


Here you're only using the first 50 columns of X. You could rewrite it as:

f_2 &lt;- numeric(n_sample)
W &lt;- 1 + X[, 1:50]/50
for (d in seq(50)) {
    f_2 &lt;- f_2 + W[,d]
}


But in the latter, you're summing along the rows of W. So you could ditch the for loop altogether:

W &lt;- 1 + X[, 1:50] / 50
f_2 &lt;- rowSums(W)


This gives us:

compute_y &lt;- function(X, A, B, epsilon) {
  n_sample &lt;- nrow(X)

  f_1 &lt;- numeric(n_sample)

  for (d in seq(100)){
    part1 &lt;- 1 + (-1)^d*A[d]*X[,d] + B[d]*sin(6*X[,d])
    f_1 &lt;- f_1 + part1
  }

  f_2 &lt;- rowSums(1 + X[, 1:50] / 50)

  f_1 * f_2 + epsilon
}


There is a way to replace the for-loop that computes f_1. 

First note you're adding 1 to f_1 one hundred times, so you might as well start with f_1 storing the value 100

f_1 &lt;- rep(100, n_sample)

for (d in seq(100)){
  part1 &lt;- (-1)^d*A[d]*X[,d] + B[d]*sin(6*X[,d])
  f_1 &lt;- f_1 + part1
}


For speed, I'll just show you how to do it:

tX &lt;- t(X)
a &lt;- colSums(c(-1, 1) * A * tX)
b &lt;- colSums(B * sin(6 * tX))
f_1 &lt;- 100 + a + b


That code would be a bit faster, but I don't think it looks as clean as your definition of f_1.

If you want you can move the code that defines X, A, B, and epsilon into a model-definition function.