Following the first direct detection of gravitational waves (GWs) from the merger of two (almost) equal mass black holes, the number of gravitational-wave detections have grown exponentially with increasing detector sensitivities. To date, tens of binary black hole and a handful of binary neutron star and neutron star - black hole binaries have been observed. These observations have been used to understand the underlying physical processes, both astrophysical and fundamental. Precise and accurate modeling of gravitational waves from such systems are paramount to the unbiased extraction of subtle effects in the gravitational-wave signal. Recent studies modeling a binary black hole ringdown signal have shown that including overtones in a ringdown waveform can model the signal closer to the merger. In this dissertation, we model a binary black hole ringdown signal including overtones, mirror modes, and subdominant modes. We show that the inclusion of mirror modes can further improve the match of the ringdown model with numerical relativity simulations. It is also shown that this more detailed model can more accurately recover the mass and spin of the final black hole. We also elucidate the role of different basis functions on the sphere, specifically, the effect of decomposing the waveform in spherical harmonics, which is the natural basis to use in a numerical simulation, versus spheroidal harmonics, which is the basis in which the radial and angular part of the signal separates. Information about the nuclear equation of state (EoS) is imprinted in the gravitational-wave signal from a binary neutron star in the form of tidal interactions between the two companion stars. Current analysis methods for extracting the tidal deformability or radii of a neutron star rely on the use of EoS-independent quasi-universal relations that relate the tidal deformabilities of the two individual stars with each other. Such quasi-relations are very useful in extracting the maximum information from a signal by reducing the dimensionality of the parameter space. However, by virtue of being quasi-relations, they contain systematic errors which become important for future observatories and, possibly, while stacking multiple current observations. We develop a methodology to mitigate these systematic errors for an unbiased and precise model selection among various equations of state. We show that unmodeled systematics can lead to the inference of the incorrect equation of state. Our method enables the use of rapid Bayesian model selection of the nuclear EoS using gravitational-wave observations. In addition to being probes of dense matter and high curvature regimes, transient gravitational-wave sources are also standard sirens. They can, therefore, act as a cosmic distance ladder that can map out distances in the Universe. When complemented by a redshift measurement from a gravitational-wave source, GWs can inform us of the evolutionary history of the Universe. Following this, GWs can provide a complimentary measurement of the Hubble constant, thus enabling the resolution of Hubble tension. We show that a sub-population of binary black hole sources, observable in current and future detector networks, can be localized to a volume in space that contains only a single galaxy on average. An electromagnetic follow-up of such sources can give a Hubble constant measurement at a precision that resolves the Hubble tension. Future observatories can probe beyond the nearby Universe and will be able to constrain other cosmological observables such as the dark matter energy density and the late-time evolution of dark energy. We contrast and compare the bounds that can be placed on various cosmological models using an electromagnetic counterpart to measure the redshift and a counterpart-less method where the redshift is measured using the tidal information between neutron stars in a binary. In the event that the independent measurement of the Hubble constant using gravitational-wave sources is in accordance with the supernova measurements, one can then use the two separate distance ladders to directly probe the electromagnetic and gravitational-wave luminosity distances. Current methods rely on restricting the modification to the luminosity distance to the gravitational sector and uses the redshift - luminosity distance relation to obtain the electromagnetic luminosity distances given a redshift measurement. We propose a method for the direct comparison of the two luminosity distances using a spatially coincident supernova measurement following a gravitational-wave event. This would be an independent and novel constrain on the variation of the two distances. We, additionally, argue that the same can be achieved with standardized kilonovae and place the first direct constraints on the variation of the two distances using the first multi-messenger observation of a binary neutron star merger GW170817. We, thereby, make the case for improved standardization modeling for kilonovae.