System Stability: Can we derive stability of a discrete system (Frequency domain, Z-transform) by applying analogous methods?
If you use the bilinear transform to map from s-plane poles and zeros to z-plane poles and zeros, the left-hand half-plane of the s-plane will be mapped to the inside of the unit circle on the z-plane. This preserves stability because on s-plane the requirement for stability is that poles are on the left-hand half-plane and on the z-plane it is that they are inside the unit circle. 

I did not examine too carefully your s and z-domain formulas, but indefinite integration is marginally stable with a pole at the origin of the s-plane, or on the unit circle of the z-plane. Implemented as a nonrecursive filter:

$$\text{out}[i] = \text{in}[i-N + 1] + \text{in}[i-N + 2] + \dots + \text{in}[i-1] + \text{in}[i],$$

integration over $N$ samples is stable. Implemented as a recursive filter: 

$$\text{out}[i] = \text{out}[i-1] + \text{in}[i] - \text{in}[i-N],$$

it is stable in fixed-point arithmetic, but can drift if implemented in floating point arithmetic due to accumulation of numerical error. Even then it can be stabilized by moving the pole slightly inside the unit circle, or by periodically resetting the filter state by calculating it using the nonrecursive filter. But sample and hold is about taking an input sample and keeping the output at that value for some time. That can be programmed directly and gives no stability problems.