Calculating an output of a system (Z- transform question)
For this specific example, using a partial fraction expansion (PFE) would be an overkill, so yo should better consider using the formula :

$$ 1 - z^{-N} = (1- z^{-1}) (1 + z^{-1} + z^{-2} + ... + z^{-N+1} $$

to simplify the Z-transform $Y(z)$ of your output as:

$$\begin{align}
Y(z) &amp;= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&amp;= \left( 1 + z^{-1} + z^{-2} + ... + z^{1-N} \right) \left( 1 + z^{-1} + z^{-2} + ... + z^{1-M} \right)\\
Y(z) &amp;= 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
Y(z) &amp;= \sum_{n=0}^{K} y[n] z^{-n} = 1 + c_1 z^{-1} + c_2 z^{-2} + ... + c_{N+M-2} z^{2-N-M} \\
\end{align}$$

Now as you can see, $Y(z)$ is a simple polynomial in $z^{-1}$. You can either explicitly multiply those two polynomials in parenthesis to find those coefficients $c_k$ to get th sample values $y[n]$ as indicated by: $$y[k] = c_k ~~~,~~~ k = 0,1,2,...,N+M-2$$

you can easily get those (coefficient $c_k$) values by the following convolution:
$$ y[n] = c[n] = a[n] \star b[n]$$ where $a[n]$ and $b[n]$ are all ones of length $N$ and $M$ respectively.

Now if you insist on a PFE method, the you could also do that as follows:
$$\begin{align}
Y(z) &amp;= H(z) X(z) = \frac{1 - z^{-N} } {1-z^{-1} } \frac{ 1 - z^{-M} }{ 1-z^{-1} } \\ \\
&amp;= \frac{1 - z^{-N} - z^{-M} + z^{-N-M} } { (1-z^{-1})^2 }\\ \\
&amp;= \frac{1} { (1-z^{-1})^2 } - \frac{z^{-N}} { (1-z^{-1})^2 } - \frac{ - z^{-M}  } { (1-z^{-1})^2 } + \frac{z^{-N-M} } { (1-z^{-1})^2 }\\ \\
\end{align}$$

Now as can be seen, I have separated the expression into four pieces each of which is a delayed version of the first one (on the left without a delay). So I will express $y[n]$ as a sum of those four delayed functions, denoting the first one as $f[n]$ we have :

$$y[n] = f[n] - f[n-N] - f[n-M] + f[n-M-N]$$ where the function $f[n]$ is the inverse Z- transform of 
$$ f[n] = \mathcal{Z}^{-1} \{ \frac{1}{(1-z^{-1})^2} \} $$

looking from a table or simply finding it yourself, one can find that 
$$f[n] = u[n] \star u[n]  = r[n]$$

where $r[n]$ is the following function :
$$r[n] = (n+1) u[n]$$ 

hence the result is simplified as:

$$ y[n] = r[n] - r[n-N] - r[n-M] + r[n-M-N] $$

instead of denoting the result in terms of ramp functions, you can also simplfity the shifted sum formula through its convolution operatos as:
$$ y[n] = (u[n] \star u[n]) -(u[n] \star u[n-N])-(u[n] \star u[n-M])+(u[n] \star u[n-N-M])$$ 
$$ y[n] = u[n] \star \left(  u[n] - u[n-N] - u[n-M] + u[n-N-M] \right)$$ 

denoting the finite length sequence on the right by $g[n]$, then the output can be found to be

$$y[n] = u[n] \star g[n] = \sum_{k=-\infty}^{n} g[k] = \sum_{k=0}^{n} g[k]  $$

The rightmost sum occurs assuming $g[n]$ is a causal sequence.