Oh, this can be made 3x faster trivially by computing x^p[i]*y^p[i] outside the inner loop.
If you are willing to lose a little bit of accuracy, you can get another 2x speedup by saving log(z) and writing
function f(x, y, z, p, q)
np = length(p)
nq = length(q
results = zeros(np * nq)
k = 1
logz = log(z)
for i in 1:np
xyp = x^p[i] * y^p[i]
for j in 1:nq
results[k] = xyp * exp(logz*q[j])
k += 1
end
end
end