  We use the recently calculated two--loop anomalous dimensions of
current-current operators, QCD and electroweak penguin operators to construct
the effective Hamiltonian for $\Delta S=1$ transitions beyond the leading
logarithmic approximation. We solve the renormalization group equations and
give the numerical values of Wilson coeff. functions. We propose a new
semi-phenomenological approach to hadronic matrix elements which incorporates
the data for $CP$-conserving $K \rightarrow \pi\pi$ amplitudes and allows to
determine the matrix elements of all $(V-A)\otimes (V-A)$ operators in any
renormalization scheme and do a renormalization group analysis of all hadronic
matrix elements $\langle Q_i(\mu) \rangle$. We compare critically our treatment
of these matrix elements with those given in the literature. We find in the NDR
scheme $\epe = (6.7 \pm 2.6)\times 10^{-4}$ in agreement with the experimental
findings of E731. We point out however that the increase of $\langle Q_6
\rangle$ by only a factor of two gives $\epe = (20.0 \pm 6.5)\times 10^{-4}$ in
agreement with the result of NA31. The dependence of $\epe$ on
$\Lambda_{\bar{MS}}$, $m_t$ and $\langle Q_{6,8} \rangle$ is presented.
