For simplicity of the equations, we simplify the problem to deal with means equal to zero.

Then the problem becomes

• X with a mean of 0 and median of -1700
• Y with a mean of 0 and median of -500.

A simple way is to sample from a distribution with three points $(X_1,Y_1)$ ,$(X_2,Y_2)$ and $(X_3,Y_3)$ with weights $w_1,w_2,w_3$.

You can find some weights that solve your problem by fiddling around. One solution is

$$w_1 = 0.25, w_2 = 0.5, w_3 = 0.25$$

$$\begin{array}{} (X_1,Y_1) &=& (400,-147.2911)\\ (X_2,Y_2) &=& (-1700,-500)\\ (X_3,Y_3) &=& (3000,1147.291) \end{array}$$

The strategy is to fix some numbers and then solve the remaining ones by a computer.

One useful trick is to set the weight of $w_2 = 0.5$ and $X_2 = -1700$ and $Y_2 = -500$ such that the median is fixed and the remaining equations are easy to solve.

In R the code would be

### function that needs to be 0
f <- function(Y_3) {
  X_3 = 3*10^3
  X_2 = -1700
  Y_2 = -500
  X_1 = -(X_3+2*X_2)
  Y_1 = -(Y_3+2*Y_2)
  X = c(X_1,X_2,X_2,X_3)
  Y = c(Y_1,Y_2,Y_2,Y_3)
  cor(X,Y)-0.97
}
f <- Vectorize(f)


### plot to visually see the maximum and minimum
Y_3 = seq(-10^4,10^4,1)
plot(Y_3, f(Y_3))

### root finding function
Y_3 = uniroot(f, c(832,10^4))$root



### these are the values
X_3 = 3*10^3
Y_3 ### this one is computed by the root finding
X_2 = -1700
Y_2 = -500
X_1 = -(X_3+2*X_2)
Y_1 = -(Y_3+2*Y_2)

The sampling is finally done by sampling the categorical distribution. With weights $w_1 = w_3 = 0.25$ and $w_2 = 0.5$