Multiscale dynamics are ubiquitous in applications of modern science. Because of time scale separation between a relatively small set of slowly evolving variables and a (typically) much larger set of rapidly changing variables, direct numerical simulations of such systems often require a relatively small time discretization step to resolve fast dynamics, which, in turn, increases computational expense. As a result, it became a popular approach in applications to develop a closed approximate model for slow variables alone, which both effectively reduces the dimension of the phase space of dynamics, as well as allows for a longer time discretization step. In this work we develop a new method for the approximate reduced model, which is based on the linear fluctuation-dissipation theorem applied to statistical states of the fast variables and designed for quadratically nonlinear and multiplicative coupling. We show that, for the two-scale Lorenz 96 model with quadratically nonlinear and multiplicative coupling in both slow and fast variables, this method produces comparable statistics to what is exhibited by an original multiscale model. In contrast, it is observed that the results from the simplified closed model with a constant coupling term parameterization are consistently less precise.