  In this paper we present a calculation of the $\Delta S=1$ effective weak
Hamiltonian including next-to-leading order QCD and QED corrections. At a scale
$\mu$ of the order of few GeV, the Wilson coefficients of the operators are
given in terms of the renormalization group evolution matrix and of the
coefficients computed at a large scale $\sim M_W$. The expression of the
evolution matrix is derived from the two-loop anomalous dimension matrix which
governs the mixing of the relevant current-current and penguin operators,
renormalized in some given regularization scheme. We have computed the
anomalous dimension matrix up to and including order $\alphas^2$ and $\alphae
\alphas$ in two different renormalization schemes, NDR and HV, with consistent
results. We give many details on the calculation of the anomalous dimension
matrix at two loops, on the determination of the Wilson coefficients at the
scale $M_W$ and of their evolution from $M_W$ to $\mu$. We also discuss the
dependence of the Wilson coefficients/operators on the regularization scheme.
