One of the main challenges facing designers of high-speed integrated circuits and interconnects is predicting the effect of the variability of geometrical and physical parameters on the circuit performance. In this thesis, an algorithm based on Parametric Model Order Reduction (PMOR) is presented for statistical analysis of large microwave and high-speed circuits with multiple stochastic parameters. Using the proposed algorithm, a set of local reduced-order parameterized circuits are derived based on adaptive frequency sampling and implicit multi-moment matching projection techniques. The local models preserve the stochastic parameters as symbolic quantities. As a result, stochastic response of the circuit can be obtained by simulating the local reduced models instead of the original large system leading to significant reduction in the computational cost compared to traditional Monte-Carlo techniques. In addition, a method is presented for high-dimensional variability analysis. This algorithm is based on two main concepts, namely node tearing for parameter partitioning and sparse grid interpolation. Node-tearing is used to localize the parameters and thus reducing the number of stochastic parameters within the subcircuits and sparse grids reduce the required number of samples for a targeted accuracy. MC analysis of the overall circuit is carried out using interface equations of much smaller dimension than the original circuit equations. On the other hand, the conventional time-domain simulation algorithms are another bottleneck for uncertainty quantification. In this thesis, a higher order time-domain algorithm based on numerical inversion of the Laplace transform (NILT) is presented. The proposed method, labeled NILTn, shows that for the same computational cost of the conventional NILT, referred to as NILT0, the approximation error is reduced by a significant factor, permitting time-marching with much larger time-steps and, consequently, saving significant computational cost per time step. In addition, while the proposed approach enables increasing the simulation time step, it highlighted the need for robust approach that can recover the full waveform in between the time points generated by NILTn. Eventually, an error estimation algorithm is presented to estimate the approximation error in a given time step.