Save computational burden by calculating the output of decimate filter at output rate?

Is it possible to write a code, in a for-loop style, that only calculates the decimated data, and outperforms the performance of built-in functions?

Yes.

(the way you asked this, here would be where the answer is over, but it's not what you meant, sooooo let's continue)

"for loop style" isn't very specific, and any algorithm that implements some form of convolution will iterate over elements, but sure:

For decimation by a factor of $M$: Decompose your filter in $M$ subfilters. Round-robin the input samples into the individual subfilters. Every $M$ input samples, calculate the output as sum of the outputs of the subfilters.

There's many publications on polyphase filters. harris' multirate signal processing is the go-to reference; but there's really many other references, too, see for example this slide deck from the Binghamton University (I'm not even sure the graphics are from that slide deck's author, my apologies if I'm wrong, but I feel I've seen the graphics before).

I recently added a different answer on the same topic, with a bit more background on the polyphase "Magic" that allows the lower rate.

The built-in filter functions provided by Matlab and a lot of realtime DSP libraries are optimized to have best computation and storage performance, but they don't provide the capability of calculating the filter output only at output sampling rate.

That's not true; any real-time DSP framework I can think of has decimating / resampling filters that do exactly that. (example: GNU Radio's decimating FIR filter function; Liquid DSP's firdecim, csdr fir_decimate, the list goes on: this is really a core functionality of every real-time signal processing toolkit.)

I don't know how efficient matlab's implementation is, but upfirdn with an interpolation ratio of 1 and a decimation > 1 should do exactly that.