Green's function retrieval by seismic interferometry (SI) exists in a variety of forms. Many of the applications of SI allow the creation of new seismic traces by crosscorrelating a wavefield recorded at two separate locations. The sum of this operation over multiple sources results in the creation of a new signal such that one of the recording locations acts as a virtual source to the other. We demonstrate that it is possible to use teleseismic shear wave transmission responses to create virtual source records for lithospheric imaging through SI. We utilize teleseismic shear wave phases for primarily two reasons: namely, (a) at large offsets the incoming wavefront approximates a plane wave due to spread of the wavefront in relation to the comparatively finite seismic array, and (b) the large distances act as a natural temporal filter to separate the incident P and S arrivals. We develop a method that allows such shear wave phases to image the subsurface. A series of forward modelling experiments with 2D elastic propagation were conducted. SP converted energy maps primarily to the vertical component. At large ray parameters, decomposition was capable of separating the P and S fields. Sufficient illumination by many ray parameters is essential to clearly resolve subsurface features; however, decimation studies indicated that it is still possible to identify strong reflectors with suboptimal sampling. Attempts to mitigate the complex effective wavelet generated by earthquake faulting behaviour were met with difficulty; deconvolution was unsuccessful due to the non-minimum phase nature of the effective source wavelet, but spectral whitening is a necessary preprocessing step to SI was necessary. Correlation panel filtering techniques, including singular value decomposition and wavenumber filtering, were tested to remove spurious energy present in the virtual source records created through SI. Spurious terms introduced by the long duration earthquake-generated wavelet were further suppressed in the common-offset domain. After creating virtual source records at every model receiver location, imaging was performed. A conventional migration technique was applied to each virtual source record and stacked to produce an image of the model subsurface. All major features of the model were resolved with a high degree of accuracy using the approach developed here. Finally, a field dataset was selected from southern Mexico with sufficient S phase sampling and prepared for imaging. By azimuthally limiting a window around the relatively linear orientation of the seismometer array, earthquakes that vary in distance correspond to sampling of unique ray parameters. A total of 43 S phase transmission responses were selected that met our criteria for magnitude and location. Geometric and topographic corrections were applied to each transmission response. The processing steps tested on our forward models were subsequently applied to this dataset and virtual source records were created via SI. The virtual source records are then migrated to form an image of the lithosphere. The resulting image of the Cocos subduction zone reproduces features present in the literature and also reveals new details at greater depths; this demonstrates the utility of the method developed within.