A fork–join system with Pareto service time distribution is considered as a model of a parallel structure of data processing. For an approximation of the mean response time of the system and its standard deviation, we apply the approach based on a combination of simulation with linear regression and the method of nonlinear Nelder–Mead optimization. Previously, no analysis of fork–join queueing systems with subsystems was carried out due to the complexity of its implementation. Nevertheless, the approach proposed here is capable of delivering approximations of various types of the correlation coefficients of the sojourn times of subtasks in a given system. The analytic expressions derived below are shown to deliver good approximations to these characteristics, as evidenced by numerical experiments. Application of the proposed approach can be extended to systems with more involved architecture, and, in particular, to systems with non-Poisson input flow and various options of distributing of service times of tasks.